# Experimental Condensed Matter Physics

32

# **Local Probing of a Superconductor's Quasiparticles and Bosonic Excitations with a Scanning Tunnelling Microscope**

**Local Probing of a Superconductor's QPs and Bosonic Excitations with an STM**

Thomas Gozlinski

**32**

Thomas Gozlinski

Thomas Gozlinski

**Local Probing of a Superconductor's Quasiparticles and Bosonic Excitations with a Scanning Tunnelling Microscope** **Experimental Condensed Matter Physics Band 32**

Herausgeber **Physikalisches Institut** Prof. Dr. David Hunger Prof. Dr. Alexey Ustinov TT-Prof. Dr. Philip Willke Prof. Dr. Wolfgang Wernsdorfer Prof. Dr. Wulf Wulfhekel

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

# **Local Probing of a Superconductor's Quasiparticles and Bosonic Excitations with a Scanning Tunnelling Microscope**

by Thomas Gozlinski

Karlsruher Institut für Technologie Physikalisches Institut

Local Probing of a Superconductor's Quasiparticles and Bosonic Excitations with a Scanning Tunnelling Microscope

Zur Erlangung des akademischen Grades eines Doktors der Naturwissenschaften von der KIT-Fakultät für Physik des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Thomas Gozlinski, M.Sc.

Tag der mündlichen Prüfung: 11. November 2022 Referent: Prof. Dr. Wulf Wulfhekel Korreferent: Prof. Dr. Jörg Schmalian

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2023 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 2191-9925 ISBN 978-3-7315-1279-0 DOI 10.5445/KSP/1000155716

# **Contents**




# **1 Introduction**

Superconductors are materials with two very unique characteristics: They are able to transmit a current without electrical resistance and can completely expel external magnetic fields from their interior (Meißner-Ochsenfeld effect). Unfortunately, all superconducting materials we know so far, only show these properties under extreme conditions, i.e. very low temperature or large pressure, which currently limits their application to large scale science projects or complex medical imaging techniques. Their potential is, however, so great and diverse that superconductivity at ambient conditions is a long-sought holy grail of condensed matter physics. Key technologies that would follow such a discovery range from lossless electrical power transmission and storage to less consumptive transportation alternatives or powerful quantum computation [1].

The history of fundamental research in superconductivity is a prime example of how scientific progress often happens in quantum leaps and in areas where we least expect it. When we think of good conductors, we think of silver, gold or platinum. The first discovered superconductor was mercury. In 1911, Kammerlingh Onnes measured a vanishing electrical resistance below a critical temperature of 4.2 K [2–4]. While mercury is a metal at room-temperature, it is a rather bad metal when compared with gold, silver or platinum. Things became only more absurd with the discovery of the first high temperature superconductor (HTS). Following the discovery of superconductivity in mercury, scientists found superconductivity in more and more metals and metallic alloys or compounds. Bernd Matthias was convinced that he had found a pattern after which materials did or did not show superconductivity which became known as Matthias' rules [5]. Loosely speaking, they state that superconductivity favours transition metals over simple metals, a high electronic density of states (DOS), high symmetry and that it does not like oxygen, magnetism or insulating phases [6]. The cuprate HTS discovered by Bednorz and Müller in 1986 [7], an antiferromagnetic insulator with copper-oxide planes, did not break one or two, but, at least to some extent, every single one of these rules. Clearly, these types of superconductors with critical temperatures of up to 138 K [8] were fundamentally different and to this day there is no complete microscopic theory that can explain the unconventional superconductivity in these systems. But the history of theory of superconductivity is also a story of great success and the currently world-record setting high-pressurized hydrides [9, 10] are a testimony to our good understanding of conventional superconductors [11]. But why is it that "bad metals" make good superconductors? In order to understand this, we have to go back to the early days of superconductivity theory.

The 1950s, roughly 40 years after superconductivity was found in mercury, turned out to be some of the most important years for the history of superconductivity, that still influence the way we understand this phenomenon today. The isotope effect [12], i.e. the dependence of the critical temperature T<sup>c</sup> and critical field H<sup>c</sup> on the atomic mass in isotopes of the same element, led Fröhlich to believe that electron-lattice interactions are of utmost importance for the superconducting state to form [13]. He could show that although the static Coulomb interaction between two electrons is repulsive, a net attractive force between two electrons can be mediated by the dynamic interaction of electrons with the crystal lattice [13, 14]. A few years later, Cooper could show that the Fermi sea of electrons is unstable against even the smallest attractive interaction, which leads electrons to form pairs, the so-called Cooper pairs [15]. A year after, in 1957, the Nobel Prize winning work of Bardeen, Cooper and Schrieffer [16] finally tied all the threads together to a coherent microscopic theory of conventional superconductivity, that could describe the relation between microscopic parameters like the binding energy of a Cooper pair and macroscopic quantities like the critical temperature. With it, the vanishing resistance, the Meißner effect and the jump in the specific heat at T<sup>c</sup> could all be explained from fundamental interactions of electrons in solids. With the extension to a strong-coupling theory [17], we also learned that there is in principle no natural upper limit for T<sup>c</sup> [18] and that a stronger coupling between electrons and lattice vibrations does not only lead to an increased normal state conductivity due to electron scattering but can also lead to higher critical temperatures.

It was during this time period, that the Norwegian physicist Giaever developed the technique of electron tunnelling spectroscopy [19], which later earned him a Nobel Prize. From there on, tunnelling experiments, tunnelling theory and the theory of superconductivity developed side by side. Today, electron tunnelling spectroscopy can be performed locally by a scanning tunnelling microscope (STM) and state-of-the-art electronics and cryogenics enable a deep look into the electronic structure of superconductors with unprecedented energy resolution. A detailed introduction into the theory of superconductivity and electron tunnelling is given in Chap. 2, followed by a description of the experimental setup and methods in Chap. 3.

Based on recent advances in our understanding of superconductivity and related phenomena,we want to go back in time and use the capability of a state-of-the-art STM to revisit superconductors that still leave room for undiscovered things. We start with two popular representatives of the cuprate family of HTSs (Chap. 4). With granular aluminium (grAl) (Chap. 5), we then turn to a granular superconductor with a T<sup>c</sup> enhancement mechanism that is not yet understood, but that was already discovered in 1966 [20]. Lastly, we want to have a look at the conventional superconductor lead (Pb) (Chap. 6 and 7).

Antiferromagnetic spin fluctuations are one of the hottest contenders for the exchange boson in unconventional superconductors. With the recent theoretical progress in the spin-fermion model [21, 22] and tunnelling experiments revealing bosonic spectra [23–25] that show remarkable resemblence to the spin resonance seen in inelastic neutron scattering (INS) [26], we slowly reach a consensus on the pairing mechanism in cuprate superconductors. J. Jandke and P. Hlobil could show in a combined experimental and theoretical study [27, 28] that the spectral density of the pairing boson is also accessible through the inelastic part of a tunnelling current. STM provides conditions under which an effective Eliashberg function, which is proportional to the spectral density of the pairing boson, can be extracted. The extraction of the effective Eliashberg function was in that case done for an iron-based superconductor with a full gap. In this work, the aim is to apply the formalism to two cuprate superconductors with a d-wave symmetry of the superconducting gap.

GrAl recently reemerged in the superconducting qubit community as a suitable superinductor in Al/AlO<sup>x</sup> based superconducting circuits [29–31] and as a template formagnetic-field resilient qubit fabrication [32, 33]. Not only are excess quasiparticles identified as a limiting loss mechanism in these circuits, transport measurements in grAl near its metal-to-insulator transition (MIT) also hint at localized magnetic moments that might be responsible for the 1/f flux noise in Al/AlO<sup>x</sup> based junctions [34, 35]. For some reason, a thorough STM study of this material is (to our knowledge) missing to this date. With our local probe of the quasiparticles on a microscopic scale, we are sensitive to spin scattering, which should manifest as a pair of Yu-Shiba-Rusinov (YSR) states below Tc. The aim in the scope of this thesis is to provide a detailed STM study of the material with special focus on the possible origins of microwave losses.

Multiband superconductivity is still a relatively young branch of superconductivity. Magnesium diboride (MgB2) with a T<sup>c</sup> of 39 K [36] is the first superconductor for which two-band superconductivity was confirmed in 2001/2002 [37, 38]. Two-band superconductivity was also recently demonstrated for lead (Pb) [39]. In MgB2, the existence of two superconducting condensates has interesting implications for the intermediate state and the dynamics of vortices that has been labeled type-1.5 superconductivity [40]. We are interested if a similar thing happens in Pb. Bulk Pb is not expected to host vortices in its intermediate state at temperatures close to its T<sup>c</sup> of 7.2 K as it is categorized as a traditional type-I superconductor [41]. However, Pb is the elemental superconductor that is closest to the Bogomol'nyi point that separates type-I and type-II, which is demonstrated in several experiments that initiate a type-I to type-II transition through a reduction of the effective coherence length [42, 43] or an increase in the magnetic penetration depth [44, 45]. Additionally, with our base measurement temperature of 25 mK, we are far below T<sup>c</sup> and in relatively unexplored territory. This is more of an explorative study that is solely motivated by scientific curiosity. We thus find it more instructing to lead the reader through the thought process behind consecutive measurements along the way and motivate them by the newly gained insights at each waypoint.

# **2 Superconductivity**

This chapter summarizes fundamental theoretical concepts of conventional and unconventional superconductivity and gives an introduction into tunnelling theory. The first section features the Ginzburg-Landau theory of superconductivity. It is followed by a section that covers important results of the microscopic theory of conventional superconductors. The sections on multiband superconductivity and Mott insulators familiarize the reader with superconductivity under unusual circumstances. The fourth section is dedicated to unconventional superconductivity with special focus on the cuprates and a spin-fluctuation mediated pairing mechanism. In the last section, tunnelling theory, which is essential for the spectroscopic measurement of superconductors, is introduced.

### **2.1 Ginzburg-Landau Theory**

A powerful macroscopic theory to describe superconductors at temperatures T ∼ Tc, i.e. close to the transition temperature, is the *Ginzburg-Landau (GL) theory*. It is a generalization of the *Landau theory* of second-order phase transitions, that enables one to study the low-temperature ordered state in the presence of spatial inhomegeneities. In the specific case of a superconductor, this inhomogeneity can be induced by electric currents or magnetic fields. The GL theory is a scalar field theory that describes the second-order phase transition through a thermodynamic parameter ψ(r), which is called the order parameter. ψ continuously changes from zero to a finite value upon deceeding the transition (or critical) temperature Tc, like in Fig. 2.1(a). The inhomegeneity of the system is treated through the inclusion of gradient terms of this scalar field. The lowest order of scalar field theory, in which these gradient terms play a role, is a ϕ 4 -theory. For T ≈ Tc, the free-energy of the system is expanded as a power series of the order parameter up to terms of fourth power in ψ. The free energy functional reads

**Figure 2.1: Ground State in GL Theory**: (a) Temperature dependence of the order parameter. (b) Free energy landscape in the normal α > 0 and spontaneously symmetry broken state α < 0. (c) Magnetic field and order parameter profile around a vortex in a type-II superconductor. The magnetic field decays exponentially on the lengthscale λ while the order parameter recovers on the lengthscale ξ.

$$F[\psi(\mathbf{r}), \mathbf{A}(\mathbf{r})] = F\_0 + \int \mathrm{d}r^3 \left[ \alpha |\psi(\mathbf{r})|^2 + \frac{\beta}{2} |\psi(\mathbf{r})|^4 + \right. $$

$$+ \frac{1}{2m^\*} \left| (-i\hbar \nabla - q\mathbf{A}(\mathbf{r})) \psi(\mathbf{r}) \right|^2 + \frac{|\mathbf{B}(\mathbf{r})|^2}{2\mu\_0} \right]. \tag{2.1}$$

F<sup>0</sup> is the free energy of the normal state above Tc, B = ∇ × A is the magnetic field and A the corresponding magnetic vector potential. For a superconductor, the charge and mass of each charge carrier is q = 2e and m<sup>∗</sup> = 2me, because the supercurrent is carried by Cooper pairs, that are made up of two electrons. The free energy landscape is shown in Fig. 2.1(b). For α > 0, the free energy functional is minimized by a single-valued ψ. For α < 0, a global U(1) symmetry is spontaneously broken and the energy landscape has more than one global minimum. In fact, if one extends the scalar field to the complex plane, i.e. ψ = |ψ|e iφ, the energy landscape forms a *mexican hat potential*. The sign change of α at T<sup>c</sup> is ensured by writing α(t) = α0(t− 1) with t = T /Tc. The minimization of Eq. (2.1) with respect to the order parameter ψ and the vector potential A yields the two *GL equations*:

$$\frac{1}{2m^\*} (-i\hbar \nabla - q\mathbf{A})^2 \psi + \alpha \psi + \beta |\psi|^2 \psi = 0 \,, \tag{2.2}$$

$$\dot{\mathbf{j}}\_s = \frac{1}{\mu\_0} \nabla \times \mathbf{B} = i \frac{q\hbar}{2m^\*} \left( [\nabla \psi^\*] \,\psi - \psi^\* \nabla \psi \right) - \frac{q^2}{m^\*} |\psi|^2 \mathbf{A}. \tag{2.3}$$

The dependence ofr was dropped for a shorternotation. In the field free,homogeneous case (A = 0, ∇ψ = 0), the first GL equation yields the ground state

$$|\psi|^2 = |\psi\_{\infty}|^2 = -\frac{\alpha(t)}{\beta}.\tag{2.4}$$

Allowing for spatial variations of the order parameter ψ = ψ<sup>∞</sup> + δψ(r), one finds from Eq. (2.2) that the deviations from the homegeneous ground state decay on a characteristic length scale

$$
\xi(t) = \sqrt{\frac{\hbar^2}{2m^\*|\alpha(t)|}} \propto \sqrt{\frac{1}{1-t}}, \quad t < 1,\tag{2.5}
$$

which is the GL coherence length. A second characteristic length scale is obtained by taking the curl of the second GL equation (2.3), which states that magnetic fields that penetrate the superconductor are screened by induced currents and are consequently exponentially damped on the length scale

$$
\lambda(t) = \sqrt{\frac{m^\*}{\mu\_0 |\psi(t)|^2 q^2}} \propto \sqrt{\frac{1}{1 - t}}, \quad t < 1,\tag{2.6}
$$

which is called the London penetration depth. Another important implication of Eq. (2.3) is that the magnetic flux, that penetrates the superconductor, is quantized. One can obtain this result by integrating the second GL equation along a closed line far away from the region where the field lines penetrate, so that the order parameter varies only in its phase ∇ψ = iψ∇φ(r) and that the net current through this loop is zero: H dsj<sup>s</sup> = 0. The single-valuedness of the order parameter then demands

$$
\oint \mathrm{d}s \nabla \varphi = \frac{q}{\hbar} \Phi = 2\pi n. \tag{2.7}
$$

The magnetic flux Φ = H dsA is hence quantized in units of the magnetic flux quantum Φ<sup>0</sup> = h/2e.

#### **2.1.1 The Ginzburg-Landau Parameter**

It turns out out, that superconductors can be classified into two types based on their response to a magnetic field, which is determined through the ratio of the two characteristic length scales

$$\kappa(t) \equiv \frac{\lambda(t)}{\xi(t)}.\tag{2.8}$$

κ is called the GL parameter. A superconductor with κ < 1/ √ 2 is classified as type-I and a superconductor with κ > 1/ √ 2 as type-II. Through the Meißner-Ochsenfeld effect, a superconductor tries to actively expel an externally applied field by its screening currents. At a certain field strength Hc, however, the associated energy that is needed to realize a superconducting ground state with Bint = 0 is simply larger than the free energy of the normal state and the system returns to the normal state

**Figure 2.2: Intermediate and Mixed Phase**: (a) Schematic domains of a type-I superconductor in the intermediate state. (b) Schematic Abrikosov lattice of a type-II superconductor in the mixed phase.

via a first-order phase transition. For a type-I superconductor of finite dimension the phase transition is known to occur locally, which leads to an intermediate phase, in which large normal and superconducting regions coexist. This scenario is depicted in Fig. 2.2(a). For a type-II superconductor, there exists an additional Shubnikov phase at Hc,<sup>1</sup> < H < Hc,2, in which the superconductor favors flux entry in the form of homogeneously distributed flux tubes carrying a single flux quantum Φ0, each. These flux tubes are known as vortices, and the lattice they form in a type-II superconductor is called the Abrikosov lattice, which is shown in Fig. 2.2(b). Below the lower critical field Hc,<sup>1</sup> the superconductor is again in the Meißner phase and above the upper critical field Hc,<sup>2</sup> it is in the normal phase. To illustrate the meaning of the coherence length and London penetration depth, the radial field and order parameter profile of a vortex in a type-II superconductor is schematically shown in Fig. 2.1(c).

One can understand the formation of these two transitionalphases of a superconductor by taking a look at the energy cost of a domain wall. For a strong type-I superconductor, ξ ≫ λ. Therefore, a single domain wall can lead to a large region in which the energy cost for expelling the magnetic field is already large but the density of the superconducting condensate |ψ| 2 is only slowly approaching its bulk value |ψ∞| 2 . The associated energy for hosting one domain wall is positive, so the system tries to minimize the surface area of domain walls. In contrast, a strong type-II superconductor with ξ ≪ λ is less effective in screening the penetrating magnetic field. But there, the energy gain due to the superconducting condensation outweighs the energy cost from screening. As a result, the associated energy for a single domain wall is negative and the system tries to maximize the surface area of domain walls. This is realized by a close-packed array of single-flux vortices.

#### **2.1.2 Flux Pinning and Flux Flow**

The flux lines inside a superconductor are a constant source of dissipation. In the presence of an electric current density j and magnetic field B, the Lorentz force density acting on the flux lines is given by [46, p.163]

$$f = j \times \mathbf{B}.\tag{2.9}$$

The movement of vortices or larger normal conducting regions through the crystal lattice equals a normal conducting current in their core, which is responsible for the dissipation. In the *flux flow* regime, the stopping force of the vortices can be described by a viscous drag coefficient [46, p.166]

$$
\eta = \frac{j\Phi}{v\_L},
\tag{2.10}
$$

where v<sup>L</sup> is the velocity of the vortex with respect to the lattice and Φ is the flux threaded through it.

In the absence of a transport current, the magnetic field lines enter a finite superconductor from the sides and are propelled to the inside by the magnetic field gradient. Due to the movement of their normal conducting cores, their motion is damped by dissipating currents and eventually stops. An array of vortices may also be moved by an electric field in their normal core induced by magnetic induction according to ∇ × E = −∂B/∂t [47]. Upon stopping the magnetic field ramp, the vortex motion is again damped by the viscous behaviour according to Eq. (2.10), i.e. the dynamics do not stop instantaneously.

The case v<sup>L</sup> = 0 for j ̸= 0 is only reached for η → ∞. This is a mechanism called *flux pinning*. So far, only a perfect, defect free superconductor was discussed. Inhomogeneities like point defects, inclusions, twin boundaries, stacking faults or dislocations are only a few crystal defects that can act as pinning centres. In general, for effective pinning, the defect should have the dimension of ∼ ξ or ∼ λ [46, p.163]. Point defects are therefore only effective in a few high temperature superconductors, in which ξ is of the same order as the lattice constant. Extended defects like dislocations, however, can also locally suppress the order parameter in elemental superconductors. From an energetic consideration, it is favorable for the flux to thread through the defect than through the superconducting bulk next to the defect, where the order parameter is already fully recovered. This local minimum of F[A] creates a pinning force on flux lines near it. In order to unpin the flux lines, the pinning force needs to be overcome by an overcompensating Lorentz force.

#### **2.1.3 Anderson-Higgs Mechanism**

The famous *Goldstone theorem*, which was born in the context of superconductivity by the works of Nambu [48], states that any spontaneous breaking of a continuous symmetry necessitates the emergence of massless (gapless) bosons, one for each generator of a broken symmetry [49, 50]. These bosons are referred to as *Nambu-Goldstone (NG) bosons*. Famous examples of NG modes in condensed matter physics are the magnons in ferromagnetic systems as a consequence of a broken rotational invariance in spin-space and acoustic phonons in crystals as a consequence of a broken rotational and translational invariance. The most original example, however, is the NG mode in a superconductor, which appears as a consequence of a broken global U(1) symmetry.

A BCS superconductor is often casually said to spontaneously break the U(1) gauge symmetry when entering the superconducting state. This is, however, not entirely correct because a gauge symmetry is not a real, physical symmetry, but merely a choice of parametrisation of the system<sup>1</sup> [52]. Therefore, it cannot be spontaneously broken [53]. What is meant, is that the global part of a U(1) phase rotation symmetry

**Figure 2.3: Anderson-Higgs Mechanism**: (a) Spontaneous symmetry breaking of U(1) leads to an extra degree of freedom, that is a massless NG mode (phase fluctuations), and a massive Higgs mode (amplitude fluctuations) of the order parameter ψ. (b) Dispersion ω(q) of the NG mode and the Higgs mode. The NG mode is "eaten up" in the gauge of the vector potential leading to massive photons that are gapped by the plasma frequency ωp. The Higgs mode can only decay into the quasiparticle continuum for |q| > 0. Both adapted from [51].

<sup>1</sup> The difference between a physical symmetry and a gauge symmetry can be understood as the difference between the rotation of an object and the rotation of the coordinate axes spanning the reference frame. While the first has a physical consequence, the second is only a description of the system in terms of different variables.

(which formally looks very much like a global gauge transformation) of the fermionic field operators is spontaneously broken. This symmetry is, however, a real, physical symmetry and tied to the number of fermions (or volume of the Fermi surface) in the system, which is not conserved anymore when electrons start to pair [54, 55]. As a consequence of the spontaneously broken symmetry, the system acquires the superconducting gap |∆| as an order parameter and an infinite number of degenerate ground states related by a phase φ, i.e. ∆ = ∆e ˆ iφ, that is necessarily universal throughout the system in its ground state. That this phase has physical implications and is not the freedom of a gauge is reflected in the DC Josephson effect. A superconductor is a special case, in which the matter field and gauge field, that is the electromagnetic vector potential A, are coupled. Meissner effect, flux quantization, Josephson effects are all phenomena based on this principle [56].

The U(1) symmetry describes the macroscopic phase of the order parameter, which has a preferred value in the ordered state. This is shown in the mexican-hat potential in Fig. 2.3(a). It turns out, that if such a symmetry can be associated with a gauge symmetry, like in the case of the electroweak interaction, then the massless NG mode is not visible, because it is "eaten up" by the gauge field and leads to massive bosons [57–61]. This is the so-called *Anderson-Higgs mechanism*. In the context of the electroweak interaction, these massive bosons are the W and Z bosons. In the context of superconductivity, the photons inside the superconductor gain a mass, which is ultimately responsible for the Meißner effect.

To see how this works, one can write down the GL functional for T < T<sup>c</sup> and insert the general ansatz ψ = (ψ<sup>∞</sup> + δψ(r))eiφ(r) , that allows for phase and amplitude fluctuations of the order parameter. After some algebraic transformations one finds [62]

$$F[\delta\psi, \mathcal{A}] = F\_0 + \int \mathrm{d}^3 r \left[ -2\alpha\delta\psi^2 + \frac{1}{2m^\*} \left| -i\hbar \nabla \delta\psi \right|^2 \right. $$

$$-\frac{\alpha}{\beta} \frac{\hbar^2}{2m^\*} \left| \nabla\varphi - \frac{q}{\hbar} \mathcal{A} \right|^2 + \frac{1}{2\mu\_0} \left| \nabla \times \mathcal{A} \right|^2 \right]. \tag{2.11}$$

This functional is invariant under the gauge transformation

$$\mathbf{A}' := \mathbf{A} - \frac{\hbar}{q} \nabla \varphi \,, \quad \psi' := \psi \mathbf{e}^{-i\varphi} . \tag{2.12}$$

By analyzing the form of Eq. (2.11), one recognizes two constant mass terms: The mass term for the amplitude fluctuations is the first term and the mass is −2α. The mass term for the field A is the third term and the mass is −αq<sup>2</sup>/(2βm<sup>∗</sup> ). As the microscopic theory shows, the mass term of the vector potential is equal to the plasma frequency ω<sup>p</sup> = p nsq <sup>2</sup>/(ε0m<sup>∗</sup>) [57] and the amplitude fluctuations are gapped by 2∆ [63], which is the characteristic energy gap of fermions in a superconductor and real part of the order parameter. The dispersion relation of these two modes is also shown in Fig. 2.3(b). It becomes obvious, that the quasiparticle gap 2∆ becomes completely empty through the Anderson-Higgs mechanism and the consequent elevation of the NG mode. The amplitude mode, which is commonly referred to as the *Higgs mode* can only decay for q ̸= 0, because the energy transfer to quasiparticles must be at least 2∆. The phase space for the decay of the Higgs mode into quasiparticles is, however, too large to observe the Higgs mode in most superconducting systems. According to a few theoretical studies, it should be possible to shrink the number of decay channels by decreasing the gap of the Higgs mode below 2∆ in superconductors with a charge-density wave [63] or disordered superconductors close to a superconductor-to-insulator transition (SIT) [64]. The phase mode is lifted in energy to the plasma frequency ωp, which is far above the gap, at least in 3D bulk systems. In superconductors an excitation of the phase mode is, therefore, equivalent to the excitation of a plasmon.

## **2.2 Microscopic Theory of Conventional Superconductivity**

This section gives a brief overview of important results of the Bardeen-Cooper-Schrieffer (BCS) theory for conventional superconductors [16] and the Migdal-Eliashberg theory for strong coupling superconductors [17, 65, 66].

#### **2.2.1 BCS Theory**

Given some attractive electron-electron interaction, the BCS formalism describes the pairing mechanism that leads to a superconducting condensate, the BCS ground state. It starts from a Fermi liquid theory of the normal state, i.e. the direct Coulomb interaction between electrons can be absorbed into single-particle excitations ("dressed" electrons) of the Fermi sea with energy ϵ<sup>k</sup> with respect to the chemical potential µ. The scattering process with the largest phase space, that destabilizes the Fermi liquid, is the one between a pair of electrons in the initial state (k ′ ↑ , −k ′ ↓) into the unprimed final state (k ↑ , −k ↓). The four-particle scattering vertex is described by the potential Vk,k′ (cf. Fig. 2.4) and the reduced Hamiltonian in second quantization reads

$$H = \sum\_{\mathbf{k}\sigma} \epsilon\_{\mathbf{k}} c\_{\mathbf{k}\sigma}^{\dagger} c\_{\mathbf{k}\sigma} + \frac{1}{N} \sum\_{\mathbf{k}k'} V\_{\mathbf{k}k'} c\_{\mathbf{k}\uparrow}^{\dagger} c\_{-\mathbf{k}\downarrow}^{\dagger} c\_{-\mathbf{k'}\downarrow} c\_{\mathbf{k'}\uparrow}.\tag{2.13}$$

Here and in the rest of this subsection, we closely follow and adopt the equations from Ref. [62]. Bardeen, Cooper and Schrieffer realized that it is useful to reformulate a mean-field approximation of this Hamiltonian in terms of expectation values for pair-annihilation and pair-creation. This leads to the complex valued BCS gap parameter ∆<sup>k</sup> defined by the BCS equation

$$
\Delta\_{\mathbf{k}} = -\frac{1}{N} \sum\_{\mathbf{k}'} V\_{\mathbf{k}\mathbf{k}'} \langle c\_{-\mathbf{k}'\downarrow} c\_{\mathbf{k}'\uparrow} \rangle = -\frac{1}{N} \sum\_{\mathbf{k}'} V\_{\mathbf{k}\mathbf{k}'} \frac{\Delta\_{\mathbf{k}'}}{2E\_{\mathbf{k}'}} \tanh \frac{E\_{\mathbf{k}'}}{2k\_B T} \tag{2.14}
$$

with the gapped quasiparticle (QP) spectrum

$$E\_{\mathbf{k}} = \sqrt{\epsilon\_{\mathbf{k}}^2 + \Delta\_{\mathbf{k}}^2}. \tag{2.15}$$

∆<sup>k</sup> describes the mean binding energy of a pair of electrons with opposite momenta and spin. Its complex conjugate describes the mean binding energy of a pair of holes with opposite momenta and spin. Bogoliubov [67] and Valatin [68] could show independently that the mean field Hamiltonian acquired by Bardeen, Cooper and Schrieffer can be diagonalized by a linear transformation that mixes electron and hole excitations. It then simplifies to

$$H = E\_{\rm BCS} + \sum\_{\mathbf{k}\sigma} E\_{\mathbf{k}} \gamma\_{\mathbf{k}\sigma}^{\dagger} \gamma\_{\mathbf{k}\sigma} \tag{2.16}$$

with the new fermionic creation and annihilation operators for the Bogoliubov particles (or Bogoliubons)

$$
\begin{pmatrix} \gamma\_{\mathbb{k}\uparrow} \\ \gamma\_{-\mathbb{k}\downarrow}^{\dagger} \end{pmatrix} = \begin{pmatrix} u\_{\mathbb{k}}^{\*} & -v\_{\mathbb{k}} \\ v\_{\mathbb{k}}^{\*} & u\_{\mathbb{k}} \end{pmatrix} \begin{pmatrix} c\_{\mathbb{k}\uparrow} \\ c\_{-\mathbb{k}\downarrow}^{\dagger} \end{pmatrix} \tag{2.17}
$$

and the superconducting ground state energy

$$E\_{\rm BCS} = \sum\_{\mathbf{k}} (\epsilon\_{\mathbf{k}} - E\_{\mathbf{k}} + \Delta\_{\mathbf{k}} \langle c\_{\mathbf{k}\uparrow}^{\dagger} c\_{-\mathbf{k}\downarrow}^{\dagger} \rangle). \tag{2.18}$$

The BCS ground state is thus a vacuum of Bogoliubons:

$$
\gamma\_{\mathbf{k}\sigma} \left| \psi\_{\rm BCS} \right\rangle = 0. \tag{2.19}
$$

From Eq. (2.16) it becomes obvious why Eq. (2.15) was called the QP spectrum with the characteristic gap ∆k. The superconducting gap is most notably seen in the DOS of the Bogoliubov QPs. We follow again Bardeen, Cooper and Schrieffer and assume a constant attractive pair interaction for electrons and holes in a thin shell around the Fermi energy of width 2ωD:

$$V\_{\mathbf{k}\mathbf{k}'} = \begin{cases} -V\_0, & \left|\epsilon\_{\mathbf{k}}\right|, \left|\epsilon\_{\mathbf{k}'}\right| < \omega\_D\\ 0, & \text{otherwise.} \end{cases} \tag{2.20}$$

ωD, the Debye frequency, is the natural cut-off frequency for an attraction mechanism mediated by phonons. Under this assumption, ∆<sup>k</sup> → ∆ and the DOS of Bogoliubov QPs is given by

$$\nu\_{\rm BCS}(\epsilon) = \begin{cases} \nu\_{\rm n}^F \frac{|\epsilon|}{\sqrt{\epsilon^2 - \Delta^2}} & \text{for } |\epsilon| > \Delta, \\ 0, & \text{for } |\epsilon| < \Delta. \end{cases} \tag{2.21}$$

Here, ν F n is the DOS at the Fermi energy in the normal state. In the BCS approximation (Eq. (2.20)), the gap equation (2.14) can be solved analytically and yields for T = 0

$$
\Delta(0) \approx 2\omega\_D \mathrm{e}^{-\frac{1}{\nu\_n^F V\_0}}.\tag{2.22}
$$

It is related to the critical temperature T<sup>c</sup> via the approximate relation

$$\frac{\Delta(0)}{k\_B T\_c} = 1.764.\tag{2.23}$$

Gorkov could derive the Ginzburg-Landau equations from the BCS equations and identified the position and temperature dependent gap ∆(r, T) as the order parameter ψ(r, T) of a BCS superconductor [69] with the same (1 − T /Tc) <sup>1</sup>/<sup>2</sup> dependence close to Tc.

#### **2.2.2 Bogoliubov-de Gennes Theory**

Often, one is interested in local inhomogeneities of the superconducting state. The Bogoliubov-de Gennes (BdG) formalism provides a microscopic description that enables us to study defects, like e.g. vortices. One starts from the BCS mean-field Hamiltonian and makes it spatially dependent. The QPs are described by the two-component spinor [62]

$$|\Psi\_{\mathbf{k}}\rangle = \begin{pmatrix} c\_{\mathbf{k}\uparrow}^{\dagger} \\ c\_{-\mathbf{k}\downarrow} \end{pmatrix} |\psi\_{\text{BCS}}\rangle \tag{2.24}$$

and by the Hamiltonian

$$H\_{\rm BdG}(k) = \begin{pmatrix} \epsilon\_k & -\Delta\_k \\ -\Delta\_k^\* & -\epsilon\_k \end{pmatrix} \tag{2.25}$$

with the eigenvalues ±E<sup>k</sup> from Eq. (2.15). The spatial dependence is obtained through a Fourier transformation of the single QP excitations with momentum k [62]:

$$\Psi(\mathbf{r}) = \frac{1}{\sqrt{N}} \sum\_{\mathbf{k}} \mathbf{e}^{ik\mathbf{r}} \left| \Psi\_{\mathbf{k}} \right\rangle \,, \tag{2.26}$$

$$\begin{split} H\_{\mathrm{BdG}}(\mathbf{r}) &= \frac{1}{N} \sum\_{\mathbf{k}} \mathbf{e}^{ik\mathbf{r}} H\_{\mathrm{BdG}}(\mathbf{k}) \\ &= \begin{pmatrix} H\_0(\mathbf{r}) & -\Delta(\mathbf{r}) \\ -\Delta^\*(\mathbf{r}) & -H\_0^\*(\mathbf{r}) \end{pmatrix} . \end{split} \tag{2.27}$$

The Schroedinger equation in real space

$$H\_{\rm BdG}(\mathbf{r})\Psi(\mathbf{r}) = E\Psi(\mathbf{r})\tag{2.28}$$

is termed the *Bogoliubov-de Gennes equation*. With the BdG equation, the local QP excitations can be studied under a local potential, e.g. a vector potential H<sup>0</sup> ∝ A(r), inside the superconductor or one can impose a local form of the gap ∆(r) in the vicinity of a defect.

#### **2.2.3 Migdal-Eliashberg Theory**

In this section, we want to discuss the mechanism of superconductivity in the presence of a strong-electron phonon interaction. One now has to take a closer look at the electron-electron scattering process via the exchange of a virtual phonon, as shown in the diagrammatic comparison of the interaction vertices for the BCS and Eliashberg theory in Fig. 2.4. Because this process involves time-dependent potentials, it is useful to do so using field operators in the interaction picture and employ the Green's function technique. Without going into much detail of the theory, we want to summarize in simple terms that knowing a system's single-particle Green's function G(r − r ′ , t − t ′ ) means knowing its single-particle Hamiltonian and thus the time evolution of physical observables. We already saw that the attractive electron-electron interaction is a dynamical process. Hence, it makes most sense to consider the Fourier-transformed Green's function in momentum and frequency space. For completeness, we want to introduce the Matsubara Green's function in the imaginary frequencies [70]

$$G\_{\sigma}(\mathbf{k}, i\omega\_{n}) = -\frac{1}{\sqrt{\beta \mathcal{V}}} \int\_{0}^{\beta} \mathrm{d}\tau \int\_{-\infty}^{\infty} \mathrm{d}\tau \langle \mathcal{T}\_{\tau} \Psi\_{\sigma}(\mathbf{r}, \tau) \Psi\_{\sigma}^{\dagger}(0, 0) \rangle \mathrm{e}^{i\omega\_{n}\tau - i\mathbf{k}\tau} \tag{2.29}$$

with β = 1/kBT, τ = it and the fermionic Matsubara frequencies ω<sup>n</sup> = (2n + 1)π/β. ⟨..⟩ is the grand canonical average, T<sup>τ</sup> the time ordering operator and Ψ<sup>σ</sup> and Ψ† <sup>σ</sup> are annihilation and creation operators for the fermion field with spin σ. The physically relevant retarded Green's function can always be obtained from the Matsubara (thermal) Green's function by analytic continuation of the imaginary frequencies onto the real axis:

$$G^{\rm R}(\mathbf{k},\omega) = G(\mathbf{k},i\omega\_n \to \omega + i0^+). \tag{2.30}$$

In the presence of strong electron-phonon coupling, the electronic excitations and lattice vibrations cannot be treated independently. Instead, the polarisation of the lattice by fluctuations in the electron density exerts a back action on the electronic system leading to electrons and phonons with renormalized dispersion relations. The dressed electron and phonon propagators can be evaluated in a Dyson series of the normal propagators (Green's functions) G0:

$$G = G\_0 + G\_0 \Sigma G = G\_0 + G\_0 \Sigma G\_0 + G\_0 \Sigma G\_0 \Sigma G\_0 + \dots \tag{2.31}$$

The electronic and bosonic complete self-energies Σ and Π due to the electron-phonon interaction are then defined by the equations

$$\left[G(\mathbf{k}, i\omega\_n)\right]^{-1} = \left[G\_0(\mathbf{k}, i\omega\_n)\right]^{-1} - \Sigma(\mathbf{k}, i\omega\_n) \,, \tag{2.32}$$

$$\left[D(\mathbf{q}, i\Omega\_n)\right]^{-1} = \left[D\_0(\mathbf{q}, i\Omega\_n)\right]^{-1} - \Pi(\mathbf{q}, i\Omega\_n). \tag{2.33}$$

The phonon Green's function D is defined analogously to Eq. (2.29) with boson field operators expanded in plain waves of momentum q and the bosonic Matsubara frequencies Ω<sup>n</sup> = 2nπ/β.

The Dyson equation (2.31) is already complete and can be described by a diagrammatic representation. But the self-energy still consists of a number of connected Feynman diagrams and in the case of electron-phonon interaction, it is useful to sort these

**Figure 2.4: Electron-Electron Interaction:** Feynman diagrams depicting the attractive electron-electron interaction through the effective four-fermion vertex in BCS theory (left) and through the virtual exchange of a phonon with wavevector q and frequency Ω in Eliashberg theory (right).

diagrams after their weight and write the Dyson equation in the self-consistent Born approximation (SCBA)

$$G = G\_0 + G\_0 \Sigma G = G\_0 + G\_0 (\Sigma^{(1)} + \Sigma^{(2)} + \dots) G. \tag{2.34}$$

This is the form in which we want to make an important approximation. Since the movement of the electrons adiabatically follows the movement of the ions (Born-Oppenheimer theorem) the interaction vertex is not renormalized but taken as the bare electron-phonon vertex. This is *Migdal's theorem* [65]. The neglected terms are of order O( p m/M), with m being the electron and M the ion mass, so only a few percent. This is diagrammatically depicted in Eq. (2.37). Eliashberg used the same argument to cut off the SCBA Dyson series for the electron and phonon propagators after the first fully self-consistent self-energy correction [17]. This is depicted in Eq. (2.35) and (2.36).

The self-energies are thus given by

$$\Sigma \approx \Sigma\_{\rm SCBA}^{(1)} = \bigvee \dots \tag{2.38}$$

$$\Pi \approx \Pi\_{\mathrm{SCBA}}^{\mathrm{(1)}} = \cdots \bigotimes \mathbf{\bullet} \tag{2.39}$$

The last thing that is needed before the effective electron-electron interaction can be analyzed is the extension of the electronic Green's function to electron-hole space (Nambu spinor representation). It is then given by the 2 × 2 matrix [70]

$$\hat{G}(\mathbf{k},\tau) = -\begin{pmatrix} \langle \mathcal{T}\_{\tau} c\_{\mathbf{k}\uparrow}(\tau) c\_{\mathbf{k}\uparrow}^{\dagger}(0) & \langle \mathcal{T}\_{\tau} c\_{\mathbf{k}\uparrow}(\tau) c\_{-\mathbf{k}\downarrow}(0) \rangle\\ \langle \mathcal{T}\_{\tau} c\_{-\mathbf{k}\downarrow}^{\dagger}(\tau) c\_{\mathbf{k}\uparrow}^{\dagger}(0) \rangle & \langle \mathcal{T}\_{\tau} c\_{-\mathbf{k}\downarrow}(\tau) c\_{-\mathbf{k}\downarrow}^{\dagger}(0) \rangle \end{pmatrix} := \begin{pmatrix} G\_{\mathbf{k}} & F\_{\mathbf{k}}\\ \overline{F}\_{\mathbf{k}} & -G\_{-\mathbf{k}} \end{pmatrix} \tag{2.40}$$

where the diagonal entries describe the normal propagators for spin-up electrons and spin-down holes and the off-diagonal entries are Gorkov's anomalous propagators. Just like in the BCS theory these anomalous operator averages describe the pairing mechanism. They naturally appear in the evaluation of the right diagram in Fig. 2.4 where two electrons scatter by exchange of a virtual phonon and in the boson self-energy. The two equations for the fermion and boson self-energies, Eq. (2.38) and (2.39), build, together with Eq. (2.35) and (2.36), a set of coupled integral equations that can be solved numerically. These are the *Eliashberg equations*. For their explicit form, which is lengthy and physically not too insightful, the reader is referred to standard literature.

An essential ingredient of the Eliashberg equations is the defnition of the *Eliashberg function* [70]

$$\alpha^2 F(\Omega) = -\frac{1}{V^2 \nu\_{\text{n}}^F} \sum\_{q=k-k'} \delta(\epsilon\_k) \delta(\epsilon\_{k'}) \frac{|g\_{q\Omega}|^2}{\pi} \mathbf{Im} D^{\mathcal{R}}(q, \Omega) \Theta(\Omega). \tag{2.41}$$

It is a momentum-integrated phonon DOS that electrons, which scatter between points of the Fermi surface, see. It is often also called the boson spectral function or bosonic glue in order to generalize it to the exchange of a different boson. In general Eliashberg theory, the gap becomes energy dependent ∆ → ∆(ϵ) and is mainly determined through the Eliashberg function. In a conventional s-wave superconductor a sharp resonance in the Eliashberg function at Ω ∗ leads to additional fine structure in the electronic DOS

$$\nu(\epsilon) = \nu\_\text{n}^F \text{Re} \frac{|\epsilon|}{\sqrt{\epsilon^2 - \Delta(\epsilon)^2}} \tag{2.42}$$

at energies ϵ = ∆<sup>0</sup> + n · Ω ∗ , as was demonstrated both theoretically [70] and in tunnelling experiments on Pb [71–74]. By defining the dimensionless electron-phonon coupling constant [18]

$$
\lambda = 2 \int \mathrm{d}\Omega \frac{\alpha^2 F(\Omega)}{\Omega},
\tag{2.43}
$$

the consequence of the coupling strength can be understood in the BCS limit λ ≪ 1. Electrons effectively gain a mass proportional to λ and the effective scattering potential V in Eq. (2.14) is substituted by λ, which nicely illustrates how a strong electron-phonon coupling leads to a stronger electron-electron attraction. Strongcoupling superconductors, like Pb, have larger ∆(0)/T<sup>c</sup> ratios than the BCS value (Eq. (2.23)).

### **2.3 Multiband Superconductivity**

The simultaneous existence of superconductivity in more than one electronic band is not surprising per se. After all, many superconducting materials have more than one electronic band crossing the Fermi level. The interplay of several superconducting bands has,however,been mostly studiedin the model two-band superconductorMgB<sup>2</sup> [36–38] and iron-pnictide superconductors [75–77]. For elemental superconductors, supposed evidence for more than one superconducting gap was found in the heat capacity of Nb, Ta and V in the 1960s [78]. Due to contradicting experimental results in the following years [79–81], however, the influence of multiband effects in these materials has remained a controversial topic. In recent years, modern theoretical methods and higher resolution in energy gap probing experiments paved the way for the discovery of two-band superconductivity in elemental Pb [39] after its prior theoretical prediction [82].

In a simple picture, the system of a two-band superconductor can be described by a model Hamiltonian that incorporates interband and intraband pairing, as well as interband scattering [83]:

$$H = \sum\_{\mathbf{k},\sigma} \epsilon\_{a,\mathbf{k}} c\_{a,\mathbf{k},\sigma}^{\dagger} c\_{a,\mathbf{k},\sigma} + \epsilon\_{b,\mathbf{k}} c\_{b,\mathbf{k},\sigma}^{\dagger} c\_{b,\mathbf{k},\sigma}$$

$$\begin{split} &+ \sum\_{\alpha,\beta,\mathbf{k}} \Delta\_{\alpha\beta} c\_{\alpha,-\mathbf{k},\uparrow}^{\dagger} c\_{\beta,\mathbf{k},\downarrow}^{\dagger} + \text{h.c.}\\ &+ \sum\_{\mathbf{k},\sigma} \Gamma c\_{a,\mathbf{k},\sigma}^{\dagger} c\_{b,\mathbf{k},\sigma} + \text{h.c.} . \end{split} \tag{2.44}$$

ϵα,<sup>k</sup> denotes the quasiparticle dispersion in band α, c † α,k,σ is the creation operator for a QP with momentum k and spin σ in the band α, ∆α,β is the superconducting gap parameter for inter- or intraband pairing and Γ describes the QP scattering amplitude between bands a and b. A large interband scattering amplitude Γ ≲ ∆ can lower the transition temperature T<sup>c</sup> in a linear manner [84], a small Γ ≪ ∆ enables one to effectively describe the multiband superconductor in a single-band model for each band. Consequently, the DOS may then be written as a superposition of contributions from each individual band α, both in the normal and superconducting state:

$$\nu\_{s(n)}(\omega) = \sum\_{\alpha} \nu\_{s(n),\alpha}(\omega). \tag{2.45}$$

In a tunnelling experiment, the individual band DOS may be probed with different spectral weight due to the momentum and energy dependence of the tunnelling matrix element. On a two-band superconductor this is notably seen as a discrepancy in the differential conductance at the quasi-particle coherence peaks [38, 39].

In order to determine how large the interband scattering amplitude is, one can look for hybridization-induced high energy gaps in the superconducting state, that have been predicted by Komendová et al. [85]. These gaps are a consequence of avoided band crossings of the hybridized bands at higher energies, hence their energetic position depends on the normal state dispersion. This exotic odd-ω pairing state can only be formed if ∆<sup>1</sup> ̸= ∆<sup>2</sup> and Γ > 0. The higher Γ, the larger the width of the hybridization-induced gaps, that appear symmetrically around the Fermi energy in a conductance spectrum. On top of that, the translational symmetry breaking of a vortex in an s-wave superconductor is supposed to lead to an odd-frequency pairing mechanism inside the vortex core [86]. Another theoretically predicted consequence of a large interband coupling is a coherence length and gap parameter matching. As was shown by Ichioka et al., for growing interband coupling, the vortex core radii ratio ξ (c) 1 /ξ(c) <sup>2</sup> → 1 and maximum gap sizes ∆2/∆<sup>1</sup> → 1 [87].

## **2.4 Mott Insulators**

It may sound paradoxical to talk about insulators when dealing with superconductivity, but for reasons that are still not perfectly understood, materials with enhanced critical temperatures very often feature a MIT at temperatures larger than Tc. This includes organic superconductors, fullerenes, bismuthates and the two classes that shall be addressed in more detail in the following: Granular superconductors and cuprates [88, p. 63]. In these classes of superconductors, the MIT is of Mott type. This means that a transition from the metallic to insulating state is caused by a charge carrier density that is so low that Thomas Fermi screening becomes very weak and electrons are localized on individual atoms due to strong on-site Coulomb repulsion. For a better understanding, it is instructive to first consider a metal crystal like copper where the conduction electrons can be thought of as a Fermi gas, so non-interacting particles, due to the effective Thomas Fermi screening. In a gedanken experiment proposed by Mott, a Mott type MIT could be initiated in such a crystal by just pulling the ions further and further apart and by doing so effectively reducing the charge carrier density through the increase in crystal volume<sup>2</sup> . The Thomas Fermi

<sup>2</sup> This is of course under the assumption that the crystal does not undergo a structural transition and does not gain a large number of defects in the process of its expansion.

**Figure 2.5: Mott MIT**: Schematic of a generic phase diagram of superconductors with a Mott MIT. As the control parameter x, that regulates the charge carrier density, falls below x<sup>M</sup> , the material undergoes a discontinuous phase transition from metal to Mott insulator. Upon decreasing T below Tc in a region of x around x<sup>M</sup> , the material undergoes a continuous phase transition to the superconducting state.

screening length rTF in the Fermi gas, which was the initial state, is becoming larger for lower charge concentrations like r 2 TF ∝ n −1/3 . As soon as rTF becomes much larger than the interatomic distance a, the Thomas Fermi screening is ineffective. The potential wells of neighbouring ions overlap and only bound states are formed. The strong electron-electron Coulomb repulsion at each atomic site renders electronic wavefunctions localized. In this insulating state, strong electron-electron interactions are dictating the materials electronic properties and as a consequence, a theoretical description in terms of single-particle states is not possible. Instead, a theoretical understanding of the ground state is usually gained in a local description of pair-wise interactions, like in the Hubbard model with the Hamiltonian:

$$H = -t\_{ij} \sum\_{\langle ij \rangle, \sigma} (c\_{i\sigma}^{\dagger} c\_{j\sigma} + c\_{j\sigma}^{\dagger} c\_{i\sigma}) + U \sum\_{i} n\_{i\uparrow} n\_{i\downarrow}.\tag{2.46}$$

Here, c (†) iσ are the annihilation (creation) operators for a fermion with spin σ at lattice site i and n = c † c are the fermionic density operators. The first term contains the tight-binding band structure via a hopping term tij between neighbouring sites ⟨ij⟩. The second term describes the on-site Coulomb repulsion of two electrons or holes via the potential U. From this Hamiltonian, it becomes obvious how a critical point can be reached when the decrease in energy due to the kinetic term is overcompensated by an increase in potential energy due to the on-site Coulomb repulsion.

In Fig. 2.5, a highly simplified, schematic phase diagram is shown for superconductors with a Mott MIT. When the control parameter x that regulates the charge carrier density falls below x<sup>M</sup> , the material undergoes a discontinuous phase transition from metal to Mott insulator, indicated by the dashed line. Upon decreasing T below T<sup>c</sup> in a region of x around x<sup>M</sup> the material undergoes a continuous phase transition to the superconducting state, indicated by the continuous line. The small cartoon depicts how an increased screening length leads to the formation of more localized states. Granular superconductors and cuprates are two prime examples in which the superconducting phase not only exists in close proximity to a Mott MIT, it also exhibits some of its highest transition temperatures in the vicinity of it [88, p. 64].

#### **2.4.1 Granular Superconductors**

Granular Superconductors are materials in which small metallic grains are separated by an insulating barrier. The control parameter that regulates the charge carrier density is the intergrain resistance, so the resistivity or thickness of the barrier. Here, the Mott insulator is characterized by a localization of electronic wavefunctions on single grains, as soon as the intergrain coupling becomes too weak. In a sense, the granular system is a mesoscopic Mott insulator [88, p. 68], because nanosized grains play the role of single atoms. A pair of metallic grains, separated by the insulating buffer, represents a normal metal tunnel junction with resistance R and capacitance C. The electrostatic energy of such a junction can only change in quantized decrements of E<sup>c</sup> = q <sup>2</sup>/2C, corresponding to the tunnelling of a single charge q. As a consequence, a tunnelling current can only flow above a threshold voltage of

$$|U| \ge U\_{\rm CB} = \frac{Q}{2C}.\tag{2.47}$$

This constraint is called the *Coulomb blockade*. It can be overcome by thermal or quantum fluctuations. Hence, in order to induce metallic behaviour, the thermal energy kBT has to be larger than E<sup>c</sup> or the time scale on which a single charge transfer occurs through quantum fluctuations τ = RC has to become sufficiently short, i.e. E<sup>c</sup> ≤ h/τ which requires R ≤ RK, where R<sup>K</sup> = h/e<sup>2</sup> is the von Klitzing constant [89, 90].

In the superconducting state, the metallic grains become superconducting and the junctions are effective Josephson junctions with Josephson coupling E<sup>J</sup> . The tunnelling charge carriers now have the charge q = 2e. Since the number of Cooper pairs on each grain N and the phase of the BCS wavefunction φ are conjugate variables, i.e. [N, φ] = i, the charge per grain Q = N(2e) and the phase cannot be measured independently. In the case of extremely decoupled grains, E<sup>c</sup> ≪ E<sup>J</sup> , the junction is in the Coulomb blockade regime and experiences large fluctuations of the superconducting phase. In the case of extremely well coupled grains, E<sup>J</sup> ≫ Ec, the junction is in the phase blockade regime, which means the superconducting phase is locked but charge fluctuations are large.

#### **2.4.2 Cuprates**

In copper-oxide based superconductors, the cuprates, the Mott criterium is evoked on the atomic scale. There, a transition from the antiferromagnetic Mott insulator to a metallic phase is initiated through the introduction of additional charge carriers via doping. A more detailed discussion of the phase diagram of the cuprates is presented in Sec. 4.1.

### **2.5 Unconventional Superconductivity**

A superconductor is classified as unconventional if the electron pairing mechanism is not mediated by electron-phonon interactions or additional symmetries besides the global U(1) symmetry are spontaneously broken for T < T<sup>c</sup> [91, p. 3]. This unconventional pairing mechanism can lead to an anisotropic gap function ∆<sup>k</sup> that is not explained by the BCS mean-field theory. In order to explain it, it is necessary to include strong electron correlations and local interactions. Especially for the technologically interesting high-temperature copper oxide superconductors, a magnetic order is competing with a superconducting phase. In these systems, the conduction and magnetic interactions are confined to the copper oxide planes and are very sensitive to doping, so the charge carrier density. Despite the absence of a complete microscopic theory that explains the unconventional electronic pairing in these systems, considerable progress has been made in the framework of a spin-fluctuation mediated pairing mechanism. This mechanism can explain the close proximity of superconductivity and magnetism as well as the resulting order parameter symmetry in the copper oxide superconductors [92].

#### **2.5.1 Order Parameter Symmetry**

The classification of superconductors into groups based on the symmetry of their pairing amplitude proves to be a powerful concept to make good physical predictions without tedious calculations. We recall that the pairing amplitude is given by

$$\Delta\_{\sigma\_1 \sigma\_2}(k) = -\frac{1}{N} \sum\_{\mathbf{k}'} V\_{\mathbf{k} \mathbf{k}' \sigma\_1 \sigma\_2} \langle c\_{-\mathbf{k}' \sigma\_1} c\_{\mathbf{k}' \sigma\_2} \rangle \tag{2.48}$$

and is determined in the self-consistent Eq. (2.14). In this eigenvalue problem of Tc, ∆<sup>σ</sup>1σ<sup>2</sup> (k) are the eigenfunctions which need to have certain symmetry properties based on their definition. The phase transition from the normal conducting to superconducting state is a second-order phase transition and therefore tightly connected with the reduction of the underlying symmetry of the system's ground state. If the full symmetry group G describes the system's normal state, then the superconducting state belongs to a subgroup H of G [56]

$$H \subset G, \quad T < T\_c. \tag{2.49}$$

An important implication of spontaneous symmetry breaking is that the Hamiltonian that describes the time-evolution of the system is still invariant under all symmetry operations of G. Only the ground state has a reduced symmetry and is not invariant under symmetry operation corresponding to the broken part of G. Starting from the simplest consideration of a Fermi liquid, in which the interaction from Eq. (2.48) is isotropic in real space and momentum space, we only need to account for symmetry of the fermionic wavefunctions, i.e. respect the Pauli principle. The full symmetry group G would look as follows:

$$G = \underbrace{SU(2)}\_{\text{spin-space}} \times \underbrace{SO(3)}\_{\text{k-space}} \times \underbrace{U(1)}\_{\text{phase}} \times \underbrace{T}\_{\text{time}-\text{reversal}} \tag{2.50}$$

First, we want to use the fact that eigenfunctions of the system are built from a basis of an irreducible representation (irrep) of the symmetry group. A representation of the group H would be

$$\mathcal{D}\_H = \mathcal{D}\_l \otimes \mathcal{D}\_s \otimes \mathcal{D}\_T \tag{2.51}$$

where D<sup>l</sup> is an irrep of momentum space and D<sup>s</sup> and D<sup>T</sup> denote representations of spin-space and the time-reversal symmetry [93]. Analogously to the hydrogen atom, we can use the spherical harmonics Ylm(k) as the irrep of the SO(3) in momentum space and since our fermions have a spin of 1/2 we can use the spin functions χsm<sup>s</sup> , which are the eigenvectors of the 2D Hermitian Pauli matrices. The eigenfunctions to the eigenvalues Tc(l, s) can then be written as [93]

$$\Delta\_{l,s}(\mathbf{k}) = \begin{cases} \sum c\_m Y\_{lm}(\mathbf{k}) \chi\_{00} & \text{Singlet}, l = \text{even} \\ \sum c\_{mm\_s} Y\_{lm}(\mathbf{k}) \chi\_{1m\_s} & \text{Triplet}, l = \text{odd} \end{cases} \tag{2.52}$$

with complex scalar coefficients c<sup>m</sup> and cm,m<sup>s</sup> . Generally, the highest eigenvalue Tc(l, s) and its corresponding eigenfunction dominate the expansion of the order parameter. However, we can already understand the order parameter composition of eigenfunctions with different orbital or spin character from this very simple example. Of special importance is the transformation behaviour of the order parameter under parity operation P. We classify the superconductors according to their parity eigenvalue P∆(k) = ∆(−k) = (−1)<sup>l</sup>∆(k) into s-wave (l=0), p-wave (l=1), d-wave (l=2) and so on.


**Table 2.1: Residual Groups**: Spin-Singlet Pairing States with Even Parity in D4<sup>h</sup> (a) and D2<sup>h</sup> (b) [94]

(a)

<sup>∗</sup>We use Mulliken symbols., †All residual symmetry groups also contain the SU(2) of spin-space.

So far we have neglected spin-orbit interaction and the crystal field. Strong spin-orbit interaction will mix eigenstates of l and s, but we can save our basis by using the total angular momentum number j and an effective pseudo-spin α [93]. Introducing the crystal field will reduce the symmetry group in k-space from a continuous rotation group SO(3) to the discrete crystallographic point group Gc. Irreps of the crystallographic point groups are tabulated in numerous review articles on this subject [93, 94]. Since our focus will be on the cuprates, we will concentrate on the orthorhombic point group D2<sup>h</sup> and compare it to the irreps of the point group D4h, the point group of a tetragonal crystal. The irreps of these two point groups for singlet-pairing (even parity) are shown in Tab. 2.1. A BCS superconductor can only be represented in the A1g, the only irrep where the residual group is identical to the point group of the normal state except for the global U(1) symmetry. No additional symmetry is broken. The basis function for A1<sup>g</sup> is 1. This means that the gap parameter can be written as a constant, which is equivalent to the conventional s-wave pairing state. In unconventional superconductors, however, an additional

**Figure 2.6: Pairing Symmetries:** 2D sketch of several (an)isotropic gap parameter functions ∆(k) inside the first Brillouin zone. The inner radial plot shows the polar dependence of ∆ for a fixed |k| (isotropic Fermi surface).

symmetry is spontaneously broken and in most cases this symmetry breaking occurs in k-space. The symmetry of the pairing state thus has a lower symmetry in momentum space than the crystal lattice. According to the basis functions of a certain irrep from Tab. 2.1, the gap function for a certain pair state ∆(k) can be expanded in harmonics of cos(ki) <sup>3</sup> where k<sup>i</sup> with i = x, y, z are the k-vectors in three orthogonal directions inside the Brillouin zone [95]. This is analagous to the atomic orbitals, where linear combinations of spherical harmonics are used to obtain real wavefunctions that are often denoted in the cartesian basis, e.g. a d<sup>z</sup> <sup>2</sup> or dx2−y<sup>2</sup> orbital. For the s-wave state (A1g) of the tetragonal crystal with point group D4h, one obtains [56, 95]

$$
\Delta\_s(\mathbf{k}) = \Delta\_s^{(0)} + \Delta\_s^{(1)}(\cos k\_x + \cos k\_y) + \Delta\_s^{(2)}\cos k\_z + \dots \tag{2.53}
$$

and for the dx2−y<sup>2</sup> -wave state (B1g)

$$
\Delta\_{d\_{x^2-y^2}}(\mathbf{k}) = \Delta\_{d\_{x^2-y^2}}^0 (\cos k\_x - \cos k\_y) + \dots \quad . \tag{2.54}
$$

Apart from these pure pairing states, it is also possible to construct mixed pair states by adding the gap parameter times i of one representation to the gap parameter of another representation, provided that both representations are 1D [95]. An example for D4<sup>h</sup> would be

$$
\Delta\_s(\mathbf{k}) + i \Delta\_{d\_{x^2-y^2}}(\mathbf{k}).\tag{2.55}
$$

It should be noted that these mixed states break time-reversal symmetry and result in two first-order phase transitions instead of one second-order phase transition [56]. Remarkable in the point group D2<sup>h</sup> is that the s- and dx2−y<sup>2</sup> -pairing state belong to the same irrep A1g. A simple superposition of the two is therefore allowed without building a time-reversal symmetry breaking mixed state or requiring two unrelated phase transitions. A more detailed discussion on the pairing symmetry in the two cuprates studied in this work follows in Sec. 4.1.2.

### **2.5.2 Macroscopic and Microscopic Description of Magnetic Response in a Metal**

In order to understand how Cooper pair formation can be mediated by spinfluctuations, it is first important to understand how the spins of itinerant electrons collectively interact with an external magnetic field. Even though the microscopic interactions may be very intricate, from a macroscopic viewpoint, the magnetic

<sup>3</sup> Using Taylor-series expansion of the cosine and sine functions make it obvious why e.g. k 2 <sup>x</sup> − k 2 <sup>y</sup> ∼ cos k<sup>x</sup> − cos k<sup>y</sup> ≈ (1 + k 2 <sup>x</sup>/2) − (1 + k 2 <sup>y</sup>/2).

properties of such a system can be described by one function: the magnetic susceptibility χ = ∂M/∂H. It is the response function of the system to an external magnetic field. Since the perturbation may be time-dependent, it is favourable to study the Fourier-transformed magnetic susceptibility, i.e. how the system reacts to an external field of a certain frequency ω and wave-vector q:

$$\chi(\mathbf{q},\omega) = \int \int \chi(\mathbf{r},t) \mathbf{e}^{-i(\mathbf{q}\mathbf{r}-\omega t)} \mathrm{d}\mathbf{r} \mathrm{d}t.\tag{2.56}$$

Since the magnetic susceptibility has to obey the principle of causality, its real and imaginary part are not independent, but related to each other through Kramers-Kronig relations. In other words, it is sufficient to either perform a diffractive scattering experiment to probe the real part of χ or an inelastic scattering experiment to probe the imaginary part responsible for damping.

Microscopically, the magnetic moment operator Mˆ (r) per unit volume is the sum over all electron magnetic moment operators mˆ <sup>i</sup>

$$
\hat{M}(\mathbf{r}) = \sum\_{i} \hat{m}\_{i} \delta(\mathbf{r} - \mathbf{r}\_{i}).\tag{2.57}
$$

Solving the time-dependent Schroedinger equation for an adiabatically turned on external field and using statistical methods, one is able to relate microscopic and macroscopic description through a magnetization (or spin) correlation function [96, p. 15]

$$\chi\_{\mu\nu}(\mathbf{q},\omega) = \frac{i}{\hbar}\mu\_0 \int\_0^\infty \langle [\hat{M}\_\mu(\mathbf{q},t), \hat{M}\_\nu(-\mathbf{q},0)] \rangle \mathbf{e}^{i\omega t} \mathbf{d}t,\tag{2.58}$$

where ⟨Aˆ⟩ = Tr(ρAˆ) is the quantum-statistical average using the density matrix ρ and [A, ˆ Bˆ] = AˆBˆ − BˆAˆ is the operator commutator. Introducing spin fluctuations by deviation from the mean value δMˆ = Mˆ −⟨Mˆ ⟩ and using the fluctuation-dissipation theorem the longitudinal dynamic structure factor S(q, ω) can be directly linked to the imaginary part of the magnetic susceptibility [96, p. 17]:

$$S\_{\mu\mu}(\mathbf{q},\omega) = \frac{\mu\_0}{\hbar} \int\_{-\infty}^{\infty} \langle \delta \hat{M}\_{\mu}(\mathbf{q},t) \delta \hat{M}\_{\mu}(-\mathbf{q},0) \rangle \mathbf{e}^{i\omega t} \mathbf{d}t$$

$$= \frac{2}{1 - e^{-\hbar\omega/k\_BT}} \text{Im}\chi\_{\mu\mu}(\mathbf{q},\omega). \tag{2.59}$$

This structure factor is directly accessible in an INS experiment. The spin correlation function and hence the magnetic susceptibility in Eq. (2.58) diverge, when the system undergoes a phase transition to a magnetically ordered state χ(q → Q) → ∞. The wave vector at which the system is most prone to do that (highest transition temperature) is called the magnetic ordering wave vector Q. Energy dissipated into the electronic system through a magnetic interaction (cf. Eq. (2.59)) can be taken up in two ways: It can produce a spin-wave, a collective excitation of the electronic spins with momentum q and energy ℏω or produce Stoner excitations. The former can be described by a bosonic excitation (a magnon) with a single-particle dispersion E(q, ω) whereas the latter form a continuum, typically in the long wave-vector and high-energy sector, where they can be thought of as local spin flips.

#### **2.5.3 Spin-Fermion Model**

Even without an external field, the itinerant electrons feel an effective magnetic field at finite temperature due to the constant temperature driven fluctuations of magnetic moments. Even though the valence electrons try to quickly screen the local spin polarizations, the process is not immediate and a retarded interaction between electrons and the background magnetic field is the consequence. Especially hot electrons, that carry enough energy and momentum to emit a magnon, can effectively attract another electron further away that feels the polarization of the environment and ultimately resolves it by absorbing the magnon. At low temperatures, this process happens only virtually (off shell) involving electrons and holes in a small energy window around the Fermi edge. The corresponding Feynman diagram to this coupling mechanism is shown in Fig. 2.7. The spin-fermion model aims to provide an effective Hamiltonian that describes the interaction between low-energy electrons and holes and collective spin-excitations. A simple Hamiltonian for this low-energy regime, assuming a short-range Hubbard-type interaction between fermions and spin fluctuations, reads [22, 97]

$$\mathcal{H} = \sum\_{\mathbf{k},\alpha} v\_{\mathbf{k}} (\mathbf{k} - k\_F) c\_{\mathbf{k},\alpha}^{\dagger} c\_{\mathbf{k},\alpha} + \sum\_{\mathbf{q}} \chi\_0^{-1}(\mathbf{q}) \mathbf{S}\_{\mathbf{q}} \mathbf{S}\_{-\mathbf{q}} + g \sum\_{\mathbf{k},\mathbf{q},\alpha,\beta} c\_{\mathbf{k}+\mathbf{q},\alpha}^{\dagger} \sigma\_{\alpha,\beta} c\_{\mathbf{k},\beta} \mathbf{S}\_{-\mathbf{q}}.\tag{2.60}$$

Here c † <sup>k</sup>,α denotes a fermion creation operator, σ<sup>i</sup> are the Pauli matrices and S<sup>q</sup> describe the magnon field. g is the coupling constant between the fermionic and bosonic degrees of freedom and χ −1 0 (q) = (ξ <sup>−</sup><sup>2</sup> + (q − Q) 2 )/χ<sup>0</sup> is the static part of the spin susceptibility. ξ is the spin correlation length and Q the magnetic ordering vector. The static part of the susceptibility has Ornstein-Zernike form, which captures the short-range interaction as its Fourier transform has the form of a Yukawa potential.

A simple separation between fermionic and bosonic degrees of freedom (which Eq. (2.60) implies), like in the case of electron-phonon interaction, is generally not possible because Migdal's theorem does not apply. Vertex corrections are suppressed by the order O(vs/v<sup>F</sup> ). The major difference compared to the phonon case is that v<sup>s</sup> now describes the spin velocity instead of the sound velocity, which is of the order of the Fermi velocity v<sup>F</sup> because the spin fluctuations stem from the fermions themselves. Consequently, not only higher order vertex corrections but also a complete renormalization of the bosonic propagator due to particle-hole bubble diagrams are necessary over a wide range of frequencies [22]. The renormalization of fermionic propagators, spin-polarization operators and vertices follows the diagrammatic series from Eq. (2.35), (2.36) and (2.37) for the phonon exchange, but now higherorder diagrams, which were discarded due to Migdal's theorem before, have to be respected.

Conceptually, one could now proceed analogously to the Migdal-Eliashberg theory and use single fermion propagators that are fully renormalized by the fermionic self-energy Σ(k, Ω) (first order and higher order diagrams in Eq. (2.35)), single boson propagators that are fully renormalized by the bosonic self-energy Π(q, ω) (first order and higher order diagrams in Eq. (2.36)) and the fully renormalized coupling constant g˜ = g + ∆g (first order and higher order diagrams in Eq. (2.37)). What makes this problem so difficult is the fact that the vertex corrections ∆g contain the full fermionic and bosonic propagators. This means that all vertex corrections and self energies have to be solved self-consistently which is a near impossible task. So it is mandatory to look for reasonable justifications to terminate the corrections. This is especially the case for a strong-coupling superconductor.

For the cuprate superconductors, there is a good justification for the neglect of many high-order terms called the hot-spot theory. One makes use of the fact, that the bare spin susceptibility in these systems peaks around the antiferromagnetic ordering vector (AFV) Q and that such a AFV connects different parts of the Fermi surface (hot spots), such that the electron coupling mechanism is also dominated by the spin-fluctuations at Q. In the normal state above T<sup>c</sup> the full spin susceptibility [22]

$$\chi(\mathbf{q},\omega) = \frac{\chi\_0 \xi^2}{1 + \xi^2 (\mathbf{q} - \mathbf{Q})^2 - \Pi(\mathbf{q},\omega)}\tag{2.61}$$

is then mainly determined by Π(Q, ω) = iω/ωsf, where ωsf is the characteristic energy scale of the boson and a temperature dependent, free parameter. At q ≈ Q, the spin-fluctuation dynamics are then governed by slow modes [98] which means Migdal's theorem can be applied in this region (no vertex corrections). Additionally,

**Figure 2.7: (Para)magnon Exchange:** Feynman diagram depicting the attractive electron-electron interaction through virtual exchange of a (para)magnon.

it turns out that close to the hot spots, higher order corrections to the fermion and boson propagator only diverge logarithmically in the coupling constant, so a renormalization is possible and the expansion can be cut off after few terms. Now the coupled integral equations can be solved self-consistently, just like in the Eliashberg theory. This is shown in detail in Ref. [21, 70]. The self-energies of the anomalous Green's function Φ(k, Ω) and normal Green's function Σ(k, Ω) can lead to a nonzero gap function ∆(k, Ω) below Tc, which is now sign changing for related hot spots ∆(k + Q) = −∆(k). The role of the Eliashberg function is now played by g 2χ ′′(ω).

The dimensionless, integrated spin spectrum is given by [28]

$$\chi''(\omega) = -3\nu\_S^0 \int \mathrm{d}^d q \mathrm{Im}\chi(\mathbf{q}, \omega)/\pi,\tag{2.62}$$

where ν 0 <sup>S</sup> denotes the normal state DOS of the superconductor. The integrated spin spectrum plays the role of the bosonic excitation spectrum in tunnelling experiments as we will see in Sec. 4.2. P. Hlobil and J. Schmalian calculated this bosonic spectrum and related functions self-consistently in the spin-fermion model (see Ref. 2.8). Fig. 2.8 shows their results. They present the imaginary part of the spin susceptibility at the AFV and the spin spectrum integrated over the 2D Fermi surface of an iron-pnictide superconductor with a full gap (s±-symmetry). In the normal state, the spin susceptibility shows the expected behaviour for overdamped spin-fluctuations. In the superconducting state the spin susceptibility is drastically renormalized: The spin spectrum at the AFV develops a spin gap for ω < ωres < 2∆. While a full spin gap is present for ω < ∆ due to the gapped fermionic degrees of freedom in the superconductor, the break-up of a Cooper pair requires an energy of ω > 2∆. In the intermediate region ∆ < ω < 2∆ a spin resonance evolves at ωres that couples to the electron-like quasiparticles, creating a peak in the DOS at ∆ + ωres. The integrated spin spectrum approaches the value of the normal state for ω ∼ 5∆. One should

**Figure 2.8: Spin-Fermion Model's Spin Resonance:** Calculated spectra in the spin-fermion model for the normal state (blue) and superconducting state (red): (a) Spin spectrum at the antiferromagnetic ordering vector (AFV) showing a spin resonance at ωres below Tc. (b) Integrated spin spectrum over the 2D Brillouin zone of a fully gapped iron-pnictide superconductor. (c) Electronic DOS. Reprinted with permission from [28].

**Figure 2.9: Tunnelling in STM and Planar Junctions**: The difference in tunnelling processes due to the tunnelling geometry are schematically illustrated for an STM geometry of tip and sample (a) and a planar junction (b). In both, the occupied (initial) local density of states (LDOS) of the top electrode and the unoccupied (final) local density of states (LDOS) of the bottom electrode are shaded in red and blue respectively. The tunnelling amplitude t(r) is largest where the electronic eigenfunctions of one lead are most strongly disturbed by the potential of the opposite lead (for the STM geometry this is the point r0).

keep in mind that the spin-fermion model fails to describe nodal parts in cuprate superconductors correctly. In addition to that, a full spin-gap is not expected for them, because of the coupling between spin-fluctuations and non-gapped electrons.

### **2.6 Quantum Tunnelling in NIS Junctions**

The BCS theory was a major success in explaining the macroscopic characteristics of superconductors, i.e. the observed jump in the specific heat, the dissipationless current and the Meißner and isotope effect. Another confirmation of one of the theory's corner stones followed three years later, when Giaver [19] realized that the slope dI/dU of the tunnelling current I in normal metal-insulator-superconductor (NIS) sandwich structures under an applied voltage U resembled the BCS DOS for quasiparticle (QP) excitations. Not only did this experiment cement the existence of the superconducting gap in the DOS of the superconductor, it also paved the way for tunnelling experiments to become one of the primary techniques for investigating the low energy DOS in not only superconductors, but conducting solids in general. Therefore, historically, the exploration of the energy spectrum of superconductors and development of tunnelling theory are closely intertwined.

### **2.6.1 Theoretical Framework of the Tunnelling Geometry**

A good theoretical understanding of Giaver's observation was developed by Bardeen [99]. Describing the insulator region as a smooth potential barrier (Wentzel–Kramers–Brillouin (WKB) approximation), he could show, that an overlap of the many-particle wavefunctions from the normal metal and the superconductor leads to tunnelling normal electrons with energy eU at a rate proportional to the superconductor's QP DOS at E<sup>F</sup> + eU. In order to keep this section concise, only formulas appropriate for the interpretation of experimental data obtained in a STM will be presented. Wavefunction and derived functions of the atomically sharp tip electrode are denoted with the index t, for the sample they carry the index s.

In contrast to planar junctions, the potential barrier separating the two leads in an STM is the vacuum. This simplifies things because the vacuum is perfectly free of states that could interfere with the tunnelling electrons. The height of the potential barrier linearly changes from the work function Φ<sup>t</sup> of the tip to Φ<sup>s</sup> of the sample. The quantum mechanical tunnelling intensity |t| <sup>2</sup> and therefore also the tunnelling current

$$|I \propto |t|^2 \propto \exp(-2\kappa d) \tag{2.63}$$

is in a WKB approximation exponentially dependent on the tip-sample distance d with κ = p 2m(Φ +¯ eU)/ℏ, where Φ = (Φ ¯ <sup>s</sup> + Φt)/2 is the averaged barrier height. An exponential dependence of the tunnelling current on the applied bias voltage only arises when the bias voltage is comparable to the work function, typically several electronvolt. For small bias voltages, eU ≪ Φ¯, this effect is negligible and κ can be considered bias independent.

The atomic sharpness of the tip electrode has two important consequences: First, Tersoff and Hamann [100] could show that for small bias voltages and a metallic tip with only s-wave states, the current through the tip-sample tunnel junction is given by

$$I \propto \int\_0^{\epsilon U} \nu\_\mathbf{s}(\mathbf{r}\_0, E\_F + \epsilon) d\epsilon. \tag{2.64}$$

The tunnelling current is, hence, directly proportional to the integrated local density of states (LDOS) νs(r, E) of the sample in the energy interval [E<sup>F</sup> , E<sup>F</sup> + eU] at the tip position r0. This is schematically illustrated in the STM geometry shown in Fig. 2.9(a). The local density of occupied (initial) states of the tip electrode are shaded in red and the local density of unoccupied (final) states of the sample are shaded in blue. The total overlap is mostly determined by the overlap at r<sup>0</sup> due to the spherical wave nature of the tip state.

Second, the contribution of inelastic tunnelling events to the total current is an order of magnitude higher than in planar tunnel junctions [73]. That is because the tip breaks the lateral translation symmetry which usually enforces the conservation of in-plane momentum of the tunnelling electron [101] in planar junctions. By removal of this condition, the phase space for a bosonic excitation grows significantly [70]. The different tunnelling paths in real space are shown next to each other in Fig. 2.9 for the STM geometry (a) and the planar junction geometry (b). In the STM geometry, tunnelling happens in a conical region below the tip, although the tunnelling amplitude t(r) quickly falls off for less vertical tunnelling paths due to the smaller integral overlap of the wavefunctions of tip and sample. In a perfect planar junction, the tunnelling paths are constrained to lead-to-lead transitions that fulfill r<sup>∥</sup> → r<sup>∥</sup> + G∥, where G<sup>∥</sup> is a in-plane lattice vector, due to the conservation of in-plane momentum. The tunnelling amplitude t(r) is roughly constant for all points along the lead edges and mainly determined by vertical tunnelling paths [101].

From Eq. (2.63) and (2.64) and Fig. 2.9(a) we get a good idea of how topographic images of the sample surface are reconstructed in the constant-current mode of an STM (see also Sec. 3.1). In order to gain a better understanding of what kind of information is contained in a bias spectroscopy of the tunnelling current, i.e. a I(U) curve, let us neglect the explicit energy dependence of the tunnelling amplitude t and drop the local dependence for now. We can keep them in the back of our minds for later. First, let us only take a look at elastic tunnelling events.

#### **2.6.2 Elastic Tunnelling**

First, let us introduce the dimensionless electronic DOS

$$
\tilde{\nu}\_{\rm s/t}(\epsilon) = \frac{\nu\_{\rm s/t}(\epsilon)}{\nu\_F^{\rm s/t}} \,, \tag{2.65}
$$

where ν<sup>F</sup> denotes the DOS at the Fermi energy. The elastic tunnelling current is given by [70, 102]

$$I^{\rm el}(U) = \frac{\sigma\_0}{e} \int\_{-\infty}^{\infty} \mathrm{d}\epsilon \,\tilde{\nu}\_{\rm s}(\epsilon) \tilde{\nu}\_{\rm t}(\epsilon - eU) \left[ S(\epsilon - eU, \epsilon) - S(\epsilon, \epsilon - eU) \right],\tag{2.66}$$

$$S(\epsilon\_1, \epsilon\_2) = n\_F(\epsilon\_1) \left(1 - n\_F(\epsilon\_2)\right). \tag{2.67}$$

The conductance constant σ<sup>0</sup> = 4πe<sup>2</sup>ν t F ν s F |t el| 2 is the purely elastic conductance for a constant sample and tip DOS and n<sup>F</sup> (ϵ) = [exp((ϵ − µ)/kBT) + 1]<sup>−</sup><sup>1</sup> are the Fermi-Dirac distributions of states with energy ϵ. Additionally, we assumed the tip DOS to be constant and the elastic tunnelling matrix element |t el| <sup>2</sup> = |⟨ψt(k ′ , ϵ)|Tˆ|ψs(k, ϵ)⟩|<sup>2</sup>

**Figure 2.10: Elastic Tunnelling Processes**: Schematic of elastic tunnelling processes between tip (left) and sample (right), where blue spheres depict electrons and white spheres holes or missing electrons. (a) Elastic electron tunnelling in an NIN junction. (b) Elastic QP tunnelling and Andreev reflection in an NIS junction. (c) Josephson tunnelling in an SIS junction.

to be energy and momentum independent. The differential conductance is then given by

$$\sigma^{\rm el}(U) = \frac{\rm dI^{\rm el}(U)}{\rm dU} = -\sigma\_0 \int\_{-\infty}^{\infty} \rm d\epsilon \tilde{\nu}\_s(\epsilon) n\_F'(\epsilon - eU). \tag{2.68}$$

For T → 0, it follows n ′ F (ϵ−eU) → −δ(ϵ−eU). Therefore, for very low temperatures, the elastic conductance is directly proportional to the density of states of the sample νs(eU) at the energy eU. Else, it is described by the convolution integral in Eq. (2.68), i.e. the DOS is "smeared out" on the scale kBT. Ramping the bias voltage U and recording the differential conductance σ = dI/dU, essentially mirrors the sample DOS around the Fermi energy.

The simple elastic tunnelling process of normal electrons between a metallic tip and sample through the vacuum barrier is schematically shown in Fig. 2.10(a). If the sample is superconducting (Fig. 2.10(b)), QPs can only tunnel directly into or out of the sample for e|U| > ∆, because of the superconducting gap in the sample's QP DOS. Only for a small tunnelling barrier, so when going continuously from the tunnelling to the point-contact regime, Andreev reflections yield finite current for e|U| < ∆. During an Andreev reflection, a hole (electron) in the tip is retro-reflected as an electron (hole) with opposite momentum and spin at the boundary to the superconductor, which leads to the annihilation (creation) of a Cooper pair in the superconducting sample. If tip and sample are superconducting (Fig. 2.10(c)), then Cooper pairs can tunnel between them through the vacuum barrier, leading to sharp spike in the tunnelling current for U = 0.

#### **2.6.3 Inelastic Tunnelling**

Although mainly consisting of elastic contributions, the total conductance also contains inelastic contributions, i.e. σ = σ el + σ inel. In STM, their weight is not insignificant. As inelastic tunnelling spectroscopy (IETS) on Pb films shows, the inelastic contributions can make up 9 % of the total conductance [73, 74]. The inelastic tunnelling current is given by [70]

$$I^{\text{inel}}(U) = \frac{\sigma\_0}{e} \frac{1}{\nu\_s^F D^2} \int\_{-\infty}^{\infty} \text{d}\epsilon \int\_0^{\infty} \text{d}\omega \times \tag{2.69}$$

$$\left(\alpha^2 F(\omega)\ddot{\nu}\_s(\epsilon)\ddot{\nu}\_t(\epsilon - \omega - eU)S(\epsilon - \omega - eU, \epsilon)n\_B(\omega)\right) \tag{2.69a}$$

$$-\alpha^2 F(\omega) \tilde{\nu}\_s(\epsilon) \tilde{\nu}\_t(\epsilon - \omega - eU) S(\epsilon, \epsilon - \omega - eU) [1 + n\_B(\omega)] \tag{2.69b}$$

$$+\alpha^2 F(\omega)\tilde{\nu}\_s(\epsilon)\tilde{\nu}\_t(\epsilon+\omega-eU)S(\epsilon+\omega-eU,\epsilon)[1+n\_B(\omega)]\tag{2.69c}$$

$$-\alpha^2 F(\omega) \ddot{\nu}\_s(\epsilon) \ddot{\nu}\_t(\epsilon + \omega - eU) S(\epsilon, \epsilon + \omega - eU) n\_B(\omega) \Big) . \tag{2.69d}$$

ω is the energy of the absorbed/emitted boson with occupation number nB(ω), the Bose-Einstein distribution function. D is a characteristic cut-off energy of the boson spectrum, or in other words, the band width of off-shell states. α is the electron-boson coupling constant and F(ω) the boson spectral density. For exclusively phononic excitations, α <sup>2</sup>F(ω) is directly proportional to the Eliashberg function and will often be used synonymously, although slight differences exist (for details see Ref. [70]).

Eq. (2.69) comprises the general inelastic tunnelling processes involving one boson of energy ω. The processes are:


Processes (2.69a) and (2.69d) are essentially absent at kBT ≪ ω because the thermal equilibrium occupation of these bosonic excitations is zero. For kBT ≪ eU only one

**Figure 2.11: Signatures of Inelastic Tunnelling**: On the left side, the inelastic tunnelling process of electrons (quasiparticles) in an NIN (NIS) junction are schematically depicted like in Fig. 2.10(a,b). The excitation of a single phonon mode with spectral density α <sup>2</sup>F(ω) (orange curve) centred at ωph leads to the total tunnelling conductance σ = σ el + σ inel and its derivative on the right. The curves for the NIN/NIS junction are shown in blue/red. Partly taken from Ref. [70]

of the remaining processes, (2.69b) or (2.69c), is present, depending on the direction of the bias voltage. For positive voltage the inelastic conductance is then given by

$$\begin{split} \sigma^{\text{inel}}(U>0) &= \frac{\text{d}I^{\text{inel}}(U>0)}{\text{d}U} \\ &\approx -\frac{\sigma\_0}{\nu\_s^F D^2} \int\_0^\infty \text{d}\omega \alpha^2 F(\omega) [1 + n\_B(\omega)] \int\_{-\infty}^\infty \text{d}\varepsilon \tilde{\nu}\_s(\epsilon) \frac{\partial S(\epsilon, \epsilon + \omega - \epsilon U)}{\partial(-eU)} \end{split} \tag{2.71}$$

$$\stackrel{T\to 0}{\rightarrow} \frac{\sigma\_0}{\nu\_\text{s}^F D^2} \int\_0^{\epsilon U} d\omega \alpha^2 F(\omega) \tilde{\nu}\_\text{s} (eU - \omega). \tag{2.72}$$

After evaluating the Bose-Einstein and Fermi-Dirac distributions for zero temperature, the inelastic conductivity is given by a convolution integral of the boson spectral density and the sample DOS over the real-valued boson frequencies ω ∈ [0, eU]. In the normal conducting state (ν˜<sup>s</sup> = 1), the second derivative of the tunnelling current is then proportional to α <sup>2</sup>F(ω). In the superconducting state, the non-constant sample DOS has to be respected.

In order to develop an intuitive understanding of Eq. (2.68) and (2.72), calculated conductance signatures for normal metal-insulator-normal metal (NIN) and NIS tunnelling, involving a single Gaussian-shaped phonon mode at ωph, are shown in Fig. 2.11. Elastic plus inelastic tunnelling via emission of a phonon of frequency ωph, with the spectral density α <sup>2</sup>F(ω) shown in orange, in the normal/superconducting sample lead to the differential conductance curve shown in blue/red. For the normal conducting sample, the differential conductance σ exhibits a step at ωph and a peak at the same energy in dσ/dU due to the inelastic tunnelling contribution. For the superconducting sample, σ shows a small (gaussian broadened) copy of the coherence peak at ωph + ∆ and a change of sign at the same energy in dσ/dU. Again, this is due to the inelastic tunnelling contribution. The increase in differential conductance for ωph < eU < ωph + ∆ and the oscillations at higher energies are due to strong electron-phonon coupling effects, which alter the DOS, so the elastic conductance, of the superconductor.

This simple example already goes to show that in the presence of strong electronboson coupling and significant inelastic tunnelling into the superconductor, a determination of the Eliashberg function from just the tunnelling spectrum becomes difficult, because α <sup>2</sup>F(ω) enters the elastic and inelastic part. For large inelastic contributions in scanning tunnelling spectroscopy (STS) measurements on conventional superconductors, the Eliashberg function can simply be determined through d 2 I/dU 2 in the normal state. For negligible inelastic contributions like in planar junction experiments, it can still be determined through McMillan and Rowell's inversion procedure [72], which fits Eliashberg theory parameters iteratively until the DOS matches the tunnelling data. For unconventional superconductors, both of the above-mentioned methods fail. For one, the spectral density of the bosonic glue undergoes drastic changes from the normal to the superconducting state and for a second, the electron-boson coupling is not of phononic nature and the use of Migdal's theorem, which the Migdal-Eliashberg theory builds up on, is not justified. In Chap. 4, we will see, how scanning tunnelling spectra on unconventional superconductors can be used to qualitatively reconstruct the boson spectral density when the contribution of inelastic contributions outweighs the fine structure of the elastic conductance.

# **3 Experimental Setup andMethods**

Scanning tunnelling microscopy is a surface science technique that grants access to the local density of electronic states at a conducting sample surface. It is therefore the ideal choice to study electronic properties in real space with local information down to the atomic scale. The resolution of this technique is highly dependent on its technical realization, e.g. vibration insulation, signal amplification, and measurement temperature. For the experiments in this work, a low temperature and ultra-clean sample surface were in all cases a prerequisite. This chapter first describes working principle and operational procedures of the experimental setup, then introduces the setups of the two used scanning tunnelling microscopes including the ultra-high vacuum (UHV) system, the cryogenic system and the decoupling from various noise sources, and lastly presents sample and probe preparation.

## **3.1 Scanning Tunnelling Microscopy and Spectroscopy**

The first prototype of a scanning tunnelling device and precursor for the scanning tunnelling microscope, we know today, was the topografiner, developed by R. D. Young and co-workers between 1965 and 1971 [103]. A metallic stylus, moveable by three-axes piezo actuators, scanned the sample surface at a constant distance recreating the surface topography. However, due to the lack of sufficient vibration

**Figure 3.1: Constant Current Mode**: In the constant current mode of the STM, a feedback-loop of the tunnelling current regulates the tip-sample distance. The topography z(x, y)|I=I<sup>t</sup> contains not only a height profile but also chemical contrast due to its sensitivity to the electron density at the surface.

insulation the tip-sample distance could not be kept small enough to be in tunnelling range but the instrument was operated in field emission range. In 1981/1982 Binnig and Rohrer presented the first version of the scanning tunnelling microscope (STM) [104], a technical improvement of the topografiner, that bettered other real space imaging techniques at the time, like field ion microscopy (FIM) and scanning electron microscopy (SEM), in terms of vertical resolution while providing similar lateral resolution. In 1986, Binnig and Rohrer received the Nobel Prize in physics for their invention - a testament to how important the delicate insulation of this sensitive instrument is.

**Working principle** In an STM, an atomically sharp metal tip is held several angstroms above the sample surface, close enough for electrons to tunnel through the vacuum barrier between sample and tip. The tunnelling current, usually pA to several nA is driven by a bias voltage, which is usually much smaller than the sample and tip work function. Due to the exponential dependence between tunnelling current and tip-sample distance (WKB-approximation), the vertical resolution of the STM can be in the pm range. These small distance adjustments are realized by a piezo scanner tube. With its four-segment electrode contact, a fine motion in all three directions is possible with good linearity and its dimensions ensure a high resonance frequency. For a coarse alignment in horizontal (x) and vertical (z) direction, the sample stage and the scanner respectively are moved by stack piezo actuators that perform slip-stick motion [105].

**Electronics** In the STM designs used in this work, the bias voltage U is applied to the sample and the tip is connected to a transimpedance amplifier that converts the tunnelling current I between tip and machine ground to FPGA readable voltages. Typical amplifications range from 10<sup>7</sup> V/A to 10<sup>9</sup> V/A. The FPGA card is the heart of the real-time controller (RC) by Nanonis. It includes an analog-to-digital converter (ADC) and digital-to-analog converter (DAC) making an easy real-time control of the STM tip possible via a LabView based computer interface. In the constant-current mode, the tip-sample distance is controlled via a feed-back loop: The user can set a set-point current It, e.g. 1 nA, and the voltages on the piezo scanner are then, in real-time, adjusted continuously such that the tunnelling current always matches the set-point current. This is schematically demonstrated in Fig. 3.1. In the constant-height mode, the tip-sample distance is kept constant (feed-back loop off) at zt.

**Topographic scan image** After piezo voltages are calibrated to a tip displacement, a rectangular area on the sample surface can be scanned in constant-current-mode by applying sawtooth voltage signals to the leads of the piezo scanner responsible for x and y fine motion. A map of the value z(x, y)|I=I<sup>t</sup> at fixed bias voltage is what is called the topographic scan image. A map of the tunnelling current I(x, y)|z=z<sup>t</sup> on the other hand is obtained in constant-height mode. Usually, a topographic scan first ensures a flat surface on which a scan in constant-height mode can then be performed safely.

**Bias spectroscopy** Bias spectroscopy is the recording of an I(U) curve at a specific tip position (x, y, z). During the voltage ramp, the feed-back loop is turned off. A (line) grid spectroscopy consists of bias spectroscopies at several points defined on a (1D) 2D grid on the surface. In scanning tunnelling spectroscopy (STS), one is often interested in the first derivative of the tunnelling current with respect to bias voltage, dI/dU. Since the tunnelling current carries noise of all frequencies (especially low frequencies can be a nuisance), the use of a lock-in amplifier to filter out most noise frequencies for the cost of a little bit of energy resolution is usually favoured over a numerical derivative.

**Lock-in amplifier** The goal of this lock-in technique is to extract the first derivative of the tunnelling current with respect to the bias voltage as noise-free as possible. First, the DC tunnelling bias U<sup>0</sup> is modulated by addition of an AC voltage Umod cos(ωt). As a result, the tunnelling current also picks up an AC component at the frequency ω which is proportional to the first derivative dI/dU|U<sup>0</sup> and the amplitude of the modulation Umod. In a second step, the tunnelling current signal is demodulated in the lock-in amplifier by multiplication with a reference signal of the frequency ω and time-averaging. This averaging cancels out the signal and noise at all frequencies that are different from ω. Only the signal and noise at ω remain. The phase of the tunnelling current and the reference signal can be aligned as follows: Out of tunnelling contact, the main contribution to the AC current at ω is due to the capacitive impedance X<sup>c</sup> of the junction. Aligning the reference phase to maximize this contribution and rotating it by another 90◦ places it in line with the AC current due to ohmic resistors. The by far largest ohmic resistance is the tunnelling resistance. Hence, with this phase alignment the demodulated signal is maximally sensitive to the tunnelling current and blind to capacitive and inductive contributions.

**Energy resolution** In order to determine the maxmimum energy resolution of the Dilution-STM at 36 mK, the dI/dU feature of a DC Josephson current, or alternatively the zeroth order Andreev bound state, was measured. Even though this feature has a functional shape governed by the dissipation of Cooper pairs through inelastic

**Figure 3.2: Energy Resolution Limit**: Differential conductance feature of a DC Josephson current (black squares). The spectrum was measured in a Pb-Pb tunnelling contact at T = 36 mK and U PK mod = 5 µV. A normal distribution was fitted to the data points (blue line).

processes including thermal photons [106], we can approximate it by a delta function centred at zero voltage<sup>1</sup> . Since the Cooper pair tunnelling signature is not broadened by a Fermi distribution in the initial and final states, an analysis of this feature lets us decouple pure distribution function effects from thermal voltage fluctuations by the junction capacitance and bias voltage modulation. We prepared a Pb tip by soft dipping a tungsten tip into a Pb(111) single crystal and measured the zero energy peak in dI/dU on the same Pb(111) crystal at small tip-sample distance. The used peak modulation voltage was U PK mod = 5 µV. The bias spectroscopy is shown in Fig. 3.2. The blue line shows the fitted function of Gaussian shape with a full width at half maximum (FWHM) of

$$\text{FWHM} = \sqrt{(2U\_{\text{mod}}^{\text{RMS}})^2 + 2\ln(2)(2\sigma\_U^{\text{el}})^2}.\tag{3.1}$$

From the fit it was found thatσ el <sup>U</sup> = (9.4±0.2) µV. U RMS mod is the root mean square (RMS) value of the AC modulation related to the peak (PK) value by U RMS mod = U PK mod/ √ 2. The temperature broadening of dI/dU features can (with minor errors) also be approximated by a Gaussian distribution. The full energy resolution (FWHM) of the instrument at milli-kelvin temperatures is finally composed of the three constituents in the following way

$$
\Delta \varepsilon = \sqrt{(3.2k\_{\rm B}T)^2 + (2eU\_{\rm mod}^{\rm RMS})^2 + ((22.1 \pm 0.5)\,\upmu\text{eV})^2}.\tag{3.2}
$$

The energy resolution of the instrument below T ∼ 80 mK and a modulation voltage of U PK mod ∼ 15 µV is mainly limited by what is thought to be voltage fluctuations due to thermal charge fluctuations in opposite sides of the tunnelling junction leads.

<sup>1</sup> The spectral function P<sup>C</sup> (U) has a FWHM ≤ 2 µV, which is significantly lower than the limitation of energy resolution due to thermal voltage noise.

**Differential conductance map** A topographic scan and a map of differential conductance dI/dU|I=I<sup>t</sup> at bias voltage U can be obtained simultaneously when the lock-in amplifier is running during the scan and the integration time per pixel is long enough. Thus, differential conductance maps are usually recorded at lower tip speeds.

**Multi-pass image** Of specialimportance in this work are so-calledmulti-pass images. Maps of differential conductance at energies eU smaller than the superconducting gap ∆ are not obtainable with a turned-on feedback-loop. In that case the tip would go into direct contact with the sample as soon as the sample exhibits a full superconducting gap because of a vanishing tunnelling current. Full grid spectroscopies are time-consuming and if one is only interested in a small number of energies but a fine spatial resolution, a multi-pass image is the best choice. First, the tip scans a row of pixels forward and backward at a bias voltage larger than ∆/e and records the profile z(x, y)I=I<sup>t</sup> along this line. Then, the tip plays back the previously recorded profile at an energy eU < ∆ and records the lock-in signal. This is also advantageous compared to the constant-height mode because areas with larger height distribution are not a problem. In order for this method to work well, high mechanical stability, low thermal drift and small piezo creep are important. Vertical tip offsets in the play-back phase can either decrease the tip-sample distance (zoff < 0) for an enhanced lock-in signal or increase it (zoff > 0) for an additional safety buffer but smaller lock-in signal.

### **3.2 Low Temperature Scanning Tunnelling Microscope**

In this work, two different setups were used: the *dilution scanning tunnelling microscope (D-STM)* and the *Joule-Thomson scanning tunnelling microscope (JT-STM)* who owe their name to their central cooling system. Both systems including UHV setup, cryostat, magnetic coil and scanning tunnelling microscope are home-built after original designs in the group of Prof. Wulf Wulfhekel [107, 108]. As the D-STM was the primary setup used for all experiments except the one described in Chap. 4, and is the more intricate one of the two, it will serve as an example to explain the UHV setup, cryostat and STM design. Differences to the JT-STM are highlighted in Sec. 3.2.2.

### **3.2.1 Dilution Scanning Tunnelling Microscope**

**UHV setup** The vacuum setup follows a classical three chamber design: The three chambers - preparation chamber, STM chamber and vacuum load-lock - are placed in an L-shape in the horizontal plane. The preparation and STM chamber are connected via a gate valve and placed in line with a shared coolable, horizontal manipulator that is mounted in the preparation chamber. The vacuum load-lock chamber is mounted to the preparation chamber at an angle of 90◦ with respect to the STM chamber connection. The load-lock has a lid, that can be opened once this chamber is vented and a small vertical manipulator in the load-lock can pass samples on into the preparation chamber. Additionally, the preparation and STM chamber possess small wobble-stick manipulators for more intricate manipulation of samples in UHV.

All chambers have a connection to the *barrel*, a cylindrical container that is pumped by a roughing pump. The barrel acts as a rough pre-vacuum for each turbomolecular pump (TMP). The load-lock, preparation and STM chamber can all be pumped with separate TMPs. The TMPs are all turned off during an STM measurement due to their vibrational noise. Instead, the preparation and STM chamber are constantly pumped by ion getter pumps which work based on the principle of physisorption and chemisorption and are thus free of mechanical noise. On top of that, the preparation chamber is equipped with a titanium sublimation pump (TSP) and the STM chamber contains a non-evaporable getter (NEG) pump, a porous reactive material with a very high surface area, ideal for the adsorption of hydrogen. Like this, base pressures below 10<sup>−</sup><sup>10</sup> mbar can be reached in the preparation and STM chamber. It should be noted that the cryostat also acts as a large cryopump. The pressures are therefore mostly in the lower 10<sup>−</sup><sup>11</sup> mbar range inside the STM chamber and even another two orders of magnitude lower inside the closed 1 K shield where the STM sits. The chamber pressures are constantly monitored by ionization gauges and ion currents in the getter pumps.

Peripheral equipment mounted on the preparation chamber includes a sputter gun (differentially pumped via load-lock), a combined instrument for low energy electron diffraction (LEED) and Auger electron spectroscopy (AES), five electron-beam evaporators for molecular beam epitaxy (MBE) and a quadrupole mass spectrometer (QMS).

**Cryostat** The wet cryostat shown in Fig. 3.3(a) features a 13 L liquid helium tank below a 23 L liquid nitrogen tank. Heat transfer by conduction is kept low by thermal insulation at their mounting point on the top flange while heat transfer by radiation is kept low by vacuum surrounding and additional radiation shields. Since radiational heat transfer from a larger temperature to a lower temperature

**Figure 3.3: D-STM cryostat and STM**: (a) Cross-sectional view of the D-STM cryostat including the STM and magnetic coil. (b) Schematic depiction of the <sup>3</sup>He/4He condensation, phase separation and gas circulation. (c) 3D rendered image (cut) of the STM body. (a,c) have been reprinted from [108] with permission of AIP Publishing.

reservoir is proportional to (T 4 <sup>&</sup>gt; − T 4 <sup>&</sup>lt;) it is optimal to use several radiation shields which are equally distributed in temperature than one shield where the temperature gradient is large. This cryostat design includes five concentric radiation shields: an intermediate radiation shield (180 K), cooled by evaporating liquid nitrogen, a shield physically connected to the liquid nitrogen tank (77 K), a second intermediate radiation shield (28 K), cooled by evaporating liquid helium, a shield physically connected to the liquid helium tank (4.2 K) and a still radiation shield (1 K).

The base temperature of 25 mK in the microscope is reached with the built-in commercial dilution unit by Bluefors Cryogenics. The coolant in such a dilution refrigerator is a <sup>3</sup>He/<sup>4</sup>He mixture. Below a critical temperature of T<sup>c</sup> = 870 mK (for a <sup>3</sup>He concentration of x<sup>c</sup> = 0.67) the mixture separates into two phases - a dilute and a concentrated phase [109]. The concentrated phase is almost pure <sup>3</sup>He, the dilute phase contains a low amount of <sup>3</sup>He and the content cannot fall below x<sup>0</sup> = 6.5 %. A dilution refrigerator makes use of this constraint to the phase equilibrium by removing <sup>3</sup>He from the dilute phase. In order to satisfy the concentration constraint in the dilute phase, <sup>3</sup>He atoms from the concentrated phase are forced to cross the phase boundary in an endothermic process [109]. A schematic picture of the dilution unit in the steady state is shown in Fig. 3.3(b). The flow direction of <sup>3</sup>He is shown as light blue arrows in the gas and as white arrows in the liquid phase. In the steady state almost all of the gas mixture is condensed inside the still and mixing chamber. The concentrated phase is lighter and thus floats atop the dilute phase from which <sup>3</sup>He is constantly removed by a TMP. Much less <sup>4</sup>He than <sup>3</sup>He is removed from the dilute phase because of its significantly lower vapor pressure. The removed <sup>3</sup>He is constantly resupplied in a closed gas cycle. In order to reach the steady state, the gas mixture has to be readily precooled. It is pressed through a small capillary with the help of a compressor and then precooled by heat exchangers in the liquid nitrogen tank and liquid helium tank (see Fig. 3.3(a)). Inside the still it is again precooled by the evaporating mixture and experiences an isoenthalpic Joule-Thomson expansion which brings the mixture below T<sup>c</sup> when entering the mixing chamber.

**STM body and mechanical noise reduction** Most parts of the STM body, shown in Fig. 3.3(c), are made from copper or gold-plated copper for quick thermalization. The only real thermal connection is established by pure silver wires to the mixer. The microscope temperature is monitored by a calibrated RuO<sup>2</sup> thermometer, mounted on the STM body. 15 × 15 mm sample plates place the sample directly under the tip. The optical access to the junction in the parked position of the microscope (elevated position above coil) allows for a quick in-situ exchange of sample or tip.

In the unparked position, the microscope hangs on three soft springs with sub-Hz resonance frequency to enhance vibrational damping in the high frequency range. For insulation from low frequency mechanical noise, the whole machine, including all chambers and cryostat, sits on pneumatic vibration isolators by Newport. Vibrational noise damping from the pumps of the <sup>3</sup>He circulation is realized by decoupling bellows and Helmholtz resonators mounted on the pump line.

**Magnetic coil and heater** Of special importance in this work is the influence of magnetic fields and temperature on superconductivity. The D-STM features a superconducting NbTi-based coil, suspended in liquid helium, with a maximum field of 6.5 T. In the unparked position the STM hangs in the axial centre of the superconducting magnet giving access to the maximum magnetic field in vertical (z) direction. The electrical power to the coil is supplied through shielded copper wire up to the liquid nitrogen tank at 78 K and through high temperature superconducting wires made from rare earth barium copper oxide (REBCO) further down in the cryostat.

The measurement temperature can be increased by using a 10 kΩ resisitive heater mounted on the STM body.

### **3.2.2 Joule-Thomson Scanning Tunnelling Microscope**

The JT-STM's cryostat lacks a dilution refrigerator and instead uses the Joule-Thomson effect to liquify pure <sup>3</sup>He which, in contrast to pure <sup>4</sup>He, is not in a suprafluid phase at this temperature. The lowest achievable temperature at the time of the experiments was 630 mK. In the parked position, the STM has a direct thermal connection to the liquid helium bath enabling a quick cooldown. In the unparked position the STM hangs in the axial centre of a 3 T split-coil magnet. The STM head itself has an upside down design compared to the one in the D-STM: The z prism and scanner approach from the bottom and the samples are facing down.

## **3.3 Tip and Sample Preparation**

### **3.3.1 Tip Preparation**

**Electrochemical Etching** The ca. 3 mm long tips were electrochemically etched from a 150 µm thick tungsten wire using a version of the DC drop-off technique [110]. The setup is shown in Fig. 3.4(a). The tungsten wire is the anode and the gold ring the cathode. An alkaline solution of NaOH (c ∼ 1 mol/L) is applied to form a meniscus inside the gold ring. Upon exceeding the activation voltage of 1.43 V, soluble tungstate ions are formed at the anode in the following electrochemical reaction [110]:

$$\text{H}(\text{s}) + 2\text{OH}^- + 2\text{H}\_2\text{O} \rightarrow \text{WO}\_4^{2-} + 3\text{H}\_2(\text{g}).\tag{3.3}$$

In reality, the required voltage for the oxidation to start is slightly higher than calculated from standard potentials. In this work, a voltage between 1.7 V and 2.0 V was frequently used. An ammeter is put in series to observe the ion current. As soon as it drops too low, the solution is reapplied to the meniscus to restore the original concentration of OH<sup>−</sup> ions. Through the dissolution of tungstate ions, the wire thins at the electrolyte/tungsten interface until it is so thin that the lower end is pulled out by gravitation. The falling part of the tungsten wire is caught by a cushion of shaving foam and used as the tip since the reaction is immediately stopped for this part after the severing. It has been reported that the upper part of the wire is blunted by the continued etching in the meniscus after separation from the lower part [111].

**Flashing** After electrochemical etching, the tips were glued inside a tip-holder by silver glue and transferred into the UHV chamber. In order to remove its 5 − 10 nm thick, insulating WO<sup>3</sup> layer [112], it is locally heated by electron bombardment in several short flashes. A positive voltage of 1 kV on the tip accelerates emitted electrons from a filament onto it. The maximum emission current between filament and tip is limited to 100 mA. The voltage at the filament is turned up quickly until the maximum heating power is reached and then immediately turned off. The maximum heating power is only present for few milliseconds. These short bombardment periods ensure a very local heating at the tip by the concentrated electric field. At the heated tip end WO<sup>3</sup> reacts with W to WO<sup>2</sup> which sublimes at 800 ◦C, far below the melting point of tungsten T<sup>m</sup> = 3410 ◦C [112]. Additional benefits of flashing are the desorption of contaminants, e.g. remains of etch solution/solvents, and/or down diffusion of such to the tip shaft.

**Shaping** In preparation of the experiments from Chap. 6 and 7, the flashed tungsten tip was additionally dipped into a Au(111) crystal in a controlled manner. Soft-dipping the tip few nm into a good metal like gold in order to pick up some gold atoms is a widely used technique to modify the shape of the tip end. The procedure is repeated until a good lateral resolution and a relatively flat tip DOS is obtained. The fact that the d-orbitals of Au are completely filled provides a tip with mostly s-wave like states around the Fermi energy and the position of the d-states with respect to the Fermi energy also makes Au the noblest of all metals [113]. The former trait provides good accordance with quantum tunnelling theory in the Tersoff-Hamann approximation and a high tunnelling conductivity but an overall lower lateral resolution at small tip-sample distances due to the large radial extension of s-wave states. The latter is important for the experiments on Pb(111), described in Chap. 6 and 7 because a pure tungsten tip has higher attractive forces towards the Pb surface atoms than a gold covered tip. A bare tungsten tip often pulls out chains of Pb atoms or even picks them up during a scan. This could lead to the destruction of a flat surface area or a superconducting tip, which should both be avoided.

**Figure 3.4: Tip and Sample Preparation**: (a) Sketch of the experimental setup for electrochemical tip etching using the DC drop-off method. (b) Schematic depiction of the sputtering and annealing process. During sputtering, an ammeter between sample and ground is installed instead of the high voltage source (Uacc) that is present during annealing by electron bombardment. (c) Cartoon of the in-situ cleaving process using a top-post. The sample and top-post are glued (1) under ambient conditions. The cleave-off (2) is done in UHV at ∼ 88 K.

#### **3.3.2 Sputtering and Annealing**

**Sputtering** A clean and smooth single crystal surface for the used metals Au(111)<sup>2</sup> and Pb(111)<sup>3</sup> was obtained by several cycles of argon sputtering and subsequent annealing. During the sputtering process, the single crystal is exposed to an argon ion beam of Uacc = 4 kV accelerating voltage at an angle of ϑ = 50◦ to the surface normal (see Fig. 3.4(b)). The argon gas<sup>4</sup> is let into the preparation chamber through a needle valve up to a pressure of pprep ∼ 2 × 10<sup>−</sup><sup>8</sup> mbar (base pressure is pprep ∼ 1 × 10<sup>−</sup><sup>10</sup> mbar). The pressures can be kept relatively low since the sputter gun is differentially pumped by the load-lock. In the sputter gun, the sputter gas is ionized by free electrons emitted from a filament running at Iem = 8 mA emission current. With these parameters, the ion current between sample and machine ground is Iion ∼ 4.5 µA. Since the raster scan of the sputter gun covers a 7 × 7 mm<sup>2</sup> area and not only the sample but also other parts of the manipulator head are connected to ground, this value is to be taken with caution when its compared to a different apparatus. It here refers to the experimental setup of the D-STM's preparation chamber. The high-energy argon ions transfer their energy and momentum to the sample's surface atoms leading to the ejection of single surface atoms or clusters.

<sup>2</sup> Single crystal Au(111) (miscut angle: ±0.1 ◦, purity: 99.999% , hat shape, diameter: 6.9 mm, thickness: 3 mm) by MaTecK GmbH

<sup>3</sup> Single crystal Pb(111) (miscut angle: ±0.1 ◦, purity: 99.999%, hat shape, diameter: 8 mm, thickness: 2 mm) by MaTecK GmbH

<sup>4</sup> Argon minican 5.0 (purity≥ 99.999 %) by Linde

As a result, the crystal surface is cleaned from surface contaminants but left with larger roughness than the initial polished crystal and some argon inclusions. With the parameters mentioned, roughly 1 ML/min is removed from the surface. Typical durations of a sputtering procedure were 10 − 20 min.

**Annealing** In a subsequent annealing step, the sample is heated to a temperature below its melting point. The annealing procedure serves two purposes: It drastically enhances surface diffusion to smoothen the surface and induces the surface segregation of impurities. These impurities can then be removed in a consecutive sputtering procedure. The sample crystal was heated using the thermal radiation of a filament placed beneath it. For annealing temperatures above 300 ◦C, the sample was additionally heated by electron bombardment (see Fig. 3.4(b)). The crystal temperature can be monitored by an infrared pyrometer<sup>5</sup> and the manipulator temperature by a thermocouple. The Au(111) crystal was annealed at Tpyro ∼ 550 ◦C. The heating protocol was 2 min 45 sec of electron bombardment at Iem = 26 mA and Uacc = 300 V. The Pb(111) crystal was annealed at up to Tpyro = 260 − 280 ◦C without high voltage for 6 min.

**Hot Sputtering** A combined process, called hot sputtering, was used for Pb(111). The simultaneous soft annealing at Tpyro < 200 ◦C and sputtering can help to keep the sputtering damage low while constantly segregating contaminants to the surface.

### **3.3.3 In-situ Cleaving**

For the cuprates used in this work, several layers were cleaved off in-situ to expose a fresh, ultra-clean surface. Due to their layered structure they can easily be cleaved in the c-plane by the top-post cleaving method shown in Fig. 3.4(c). The sample is with one side glued to the sample plate and with the other side to a metal post on top of it. The adhesive is a conducting epoxy glue<sup>6</sup> . As a result, the crystal is mechanically clamped (and therefore strained after a cooldown) and the tunnelling current can flow through the conductive glue which ensures a good electrical contact. The samples were cleaved off on the precooling stage in the STM chamber of the respective machine at ∼ 88 K. The top-post was kicked off using a wobble-stick manipulator.

<sup>5</sup> OPTCTL3MH2 by Optris

<sup>6</sup> EPO-TEK H20E by EPOTEK

# **4 Bosonic Excitation Spectrum in** Bi2Sr2CaCu2O8+<sup>δ</sup> **and**YBa2Cu3O6+x

From a technological viewpoint, high temperature superconductors (HTS) represent the most desirable class of unconventional superconductors. Among them, the cuprates, which are layered copper oxides with perovskite structure, are perhaps the most famous representatives. The discovery of the first HTS lanthanum barium copper oxide (LBCO) [7] was awarded with a Nobel Prize in 1987 and soon after, other cuprate superconductors with transition temperatures above liquid nitrogen temperature followed. Prime examples are the widely studied systems of yttrium barium copper oxide (YBCO) [114] and bismuth strontium calcium copper oxide (BSCCO) [115]. In this chapter, we address the open question of what the bosonic glue in the cuprates is. After introducing the studied materials Bi2Sr2CaCu2O8+<sup>δ</sup> (Bi2212) and YBa2Cu3O6+x (Y123) and giving a brief overview of the scientific context, inelastic scanning tunnelling data is tested to tunnelling theory based on the spin-fermion model of spin-fluctuation mediated superconductivity.

## **4.1 Structural and Electronic Properties**

#### **4.1.1 Structural Properties of Y123 and Bi2212**

The high-temperature superconducting phase of Bi2212 crystallizes in an orthorhombic Amaa structure (point group D2h). The unit cell is displayed in Fig. 4.1(a) with lattice constants a = 5.41 Å, b = 5.40 Å and c = 30.72 Å. This orthorhombic unit cell contains each of its elements four times its stoichiometric value. The slight difference in a and b supposedly results from the formation of a superlattice in the Bi-O planes [118]. Along the c-axis CuO<sup>2</sup> planes are sandwiched between SrO and Ca planes. The Ca planes are mirror planes of the unit cell. The SrO-CuO2-Ca-CuO2-SrO stacks are then again sandwiched by BiO planes. The BiO planes provide only weak bonding to each other. Careful micromechanical cleaving of this crystal preferably leads to a separation between the BiO planes due to their weak van-der-Waals bonds. The

**Figure 4.1: Structure Models:** Crystal structures of (a) Bi2Sr2CaCu2O8+<sup>δ</sup> in the orthorhombic phase (space group Amaa) (b) YBa2Cu3O6+x in the orthorhombic phase (space group Pmmm). Models were created using crystallographic data from Ref. [116, 117].

appearance of other termination layers along the c-axis is possible but occurs in rarer cases [119].

The high-temperature superconducting phase ofY123 crystallizes in an orthorhombic Pmmm structure (point group D2h). The unit cell is displayed in Fig. 4.1(b) with lattice constants a = 3.82 Å, b = 3.89 Å and c = 11.68 Å. Starting from the perovskite structure of YBa2Cu3O9, one oxygen atom is removed in every third CuO<sup>2</sup> plane and one oxygen atom is removed in the Y-plane. The result can be viewed as a layered unit cell that contains two CuO<sup>2</sup> planes separated by a Y-plane and one plane of CuO chains along the b-axis which are separated from the CuO<sup>2</sup> planes by BaO planes. The composite with low oxygen concentration x < 0.35 acquires tetragonal symmetry without CuO chains. The structural transition between the tetragonal and orthorhombic phase takes place inside the range x = 0.35 − 0.65 [118]. Unlike Bi2212, Y123 has no natural cleaving plane. The experimentally found predominant termination layers after cleaving the crystal are the BaO and CuO chain layers [120].

#### **4.1.2 Unconventional Pairing in Y123 and Bi2212**

**Figure 4.2: Electronic Properties:** (a) Schematic phase diagram of hole-doped cuprate superconductors. Adapted from [22]. (b) Typical 2D Fermi surface contour of the hyperbolic bands in the copper oxide plane of cuprate superconductors including the AFV Q that connects hot spots on the Fermi surface.

Undoped cuprates tend to be antiferromagnetic insulators at low temperatures. In Bi2212 and Y123, sufficient hole doping, by regulating the oxygen concentration during growth, can destroy the magnetic order through the injection of charge carriers and reveal a rich variety of quantum matter phases. In Fig. 4.2(a), a simplified phase diagram for the hole-doped superconductors in the temperature-doping plane is shown. A superconducting dome with parabolic Tc(p) dependence, where p denotes the number of holes per Cu atom, is found to be quite universal for cuprate superconductors [122–124]. The oxygen concentration for which T<sup>c</sup> is maximal is called optimal doping, lower (higher) doped samples are considered underdoped (overdoped). For the hole doped double-layer HTS like Y123 and Bi2212, the dome is asymmetric with a steep decline of T<sup>c</sup> in the overdoped regime when the hole doping is assumed to correspond to chemical composition (p ≈ x). A possible explanation for this trend is a redistribution of oxygen in the overdoped regime [124]. The parabolic dome is, however, recovered if the hole doping is inferred from the Fermi surface volume in angle-resolved photoemission spectroscopy (ARPES) measurements [125]. For Bi2212, the dome is centred around p = 0.18 with T max <sup>c</sup> ≈ 91 K [125], for Y123 T max <sup>c</sup> ≈ 93 K is reached at p ≈ 0.17 [126]. The most mysterious and least understood part of the phase diagram is the anomalous normal state (strange metal) at high temperatures that exhibits non-Fermi liquid behaviour<sup>1</sup> . This behaviour is

<sup>1</sup> A trademark of this anomalous normal state is a resistivity that scales linearly instead of quadratically (Fermi-liquid) with temperature.

**Figure 4.3: Spin Excitation in INS:** Inelastic neutron scattering data of optimally doped Bi2212 (a-b) and Y123 (c-d): (a,c) Q scans at the indicated temperatures revealing the resonance mode at the AFV (π,π,l). (b,d) Energy scans at the AFV showing the dispersion of the spin excitation at ωres ≈ 43 meV for Bi2212 and ωres ≈ 41 meV for Y123. For low excitation energies both materials show a spin gap. (a,b) are reprinted by permission from Springer Nature Customer Service GmbH: Springer Nature [26], Copyright (1999). (c,d) are reprinted from [121], Copyright (1991), with permission from Elsevier.

Constant energy scans at excitation energy 43meV of the neutron

most pronounced in the underdoped cuprates at elevated temperatures between the superconducting and magnetically ordered phase, where the material enters a so-called pseudogap phase below the temperature T ∗ . The name is derived from the opening of a gap in the electronic spectrum that can be seen in ARPES or STM measurements. If this phase is closely tied to spin fluctuations and represents preformed Cooper pairs [127] or if it is the result of a competing order unrelated to superconductivity, like nematicity [126], remains to be clarified.

The relatively large c/a ratio in the cuprates leads to an effective two-dimensional electronic bandstructure lying in the 2D Brillouin zone (BZ) of the copper-oxide planes which are separated by ionic spacer regions. As a result, several electronic properties of the normal, but also the superconducting state, like conductivity, superconducting coherence length or magnetic penetration depth have large out-ofplane anisotropy [56]. The strong confinement of charge carriers to the copper-oxide planes gives rise to the strong doping dependence of electronic properties and the lowered dimensionality favours quantum-criticality [22, 28]. Fig. 4.2(b) shows the generic Fermi surface contour (orange line) stemming from the square CuO<sup>2</sup> plane in a cuprate superconductor in the 2D BZ. The shaded area around the Γ point represents the occupied states and the white area round the M points the unoccupied states. At zero doping (δ = 0) the percentage area of the two is identical. In the case of hole doping (δ > 0) the relative area of unoccupied states in the first BZ dominates.

Nuclear magnetic resonance (NMR) experiments in the early 1990s could confirm a spin-singlet state in the cuprates [128] and thus a superconducting state with an even parity orbital wavefunction. As was already hinted at in Sec. 2.5.1, the famous representatives of the cuprate family Y123 and Bi2212 are in an orthorhombic D2<sup>h</sup> phase at temperatures close to T<sup>c</sup> and the relevant doping rage. Hence, for the irrep A1g, a superposition of extended s-wave (s + s<sup>±</sup> + ..) and dx2−y<sup>2</sup> basis functions is technically allowed. A pairing symmetry based on the other irreps is unlikely considering the nearest-neighbour and next-nearest neighbour positions on the rectangular CuO<sup>2</sup> lattice to which the conduction is restricted. Due to the small orthorhombicity in these systems, the superconductivity is widely believed to arise from an unconventional spin-fluctuation pairing mechanism of predominantly dx2−y<sup>2</sup> character in a near tetragonal D4<sup>h</sup> phase<sup>2</sup> . There, such a dx2−y<sup>2</sup> pairing symmetry belongs to the B1<sup>g</sup> irrep and consequently lowers the total symmetry in k-space. This order parameter symmetry is readily supported by numerous experimental evidence including results from thermodynamic measurements, ARPES, NMR, Raman scattering, SQUID interferometry and many more (see Ref. [56, 92] and references therein). An extensive summary of convincing high resolution ARPES data has been published by A. Damascelli and co-workers [129]. Still, the argument that an A1<sup>g</sup> order parameter representation reveals no additional broken symmetry remains strong and consequently an s + d pairing symmetry and even an involved conventional pairing mechanism are not completely out of the picture yet [56].

The immediate vicinity of an antiferromagnetic phase close to superconductivity hints towards strong magnetic fluctuations and strong short ranged spin correlations, i.e. an enhanced magnetic susceptibility (cf. Eq. (2.58)). On the (almost) square copper oxide planes, antiferromagnetic stripe order is realized for the AFV Q ≈ (π, π). As shown in Fig. 4.2(b), this wave vector connects hot spots on the Fermi surface at

<sup>2</sup> This means the 2D symmetry group is C4<sup>v</sup> in the copper-oxide plane instead of C2<sup>v</sup> and so the CuO<sup>2</sup> unit cell is square instead of a rectangular one.

which the DOS is large. Even outside the magnetically ordered state, where the spin correlation length is smaller, the susceptibility has been shown to be maximal around Q. In fact, INS experiments show evidence of a magnetic resonance mode at the AFV and opening of a spin-gap below T<sup>c</sup> in nodal directions, as shown in Fig. 4.3<sup>3</sup> . The magnetic resonance is a collective spin-1 excitation mode in the superconducting state and is arising from the renormalization of spin-fluctuations through the gapped fermionic degrees of freedom, as could be demonstrated in the spin-fermion model (see Sec. 2.5.3). Experimentally, the resonance mode was found to be ubiquitous in the cuprate superconductor family and its energy is tied to T<sup>c</sup> and ∆ [130].

A spin-fluctuation mediated coupling mechanism in the spin-fermion model framework, as described in Sec. 2.5.3, logically combines the dx2−y<sup>2</sup> gap parameter with the hyperbolic 2D bands and magnetic susceptibility found in these materials. Since this theory is rather complicated, a simple and short answer to the question how dx2−y<sup>2</sup> pair symmetry and the AFV Q = (π/a, π/a) are related could sound as follows: Suppose an electron with momentum q hops from Cu lattice site i + 1 to an adjacent Cu lattice site i on the square CuO<sup>2</sup> lattice. Then, it essentially carries the relative phase qa which is equal to π for q = Q. Now, the on-site Coulomb repulsion on a lattice site i drastically reduces the correlation for two electrons on the same site. It can, however, become smaller for electrons with opposite signed wavefunctions that stem from an l > 0 orbital because, there, these electrons are spatially separated. The electrons from the Cu dx2−y<sup>2</sup> orbitals provide this trait and hence allow for a non-zero correlation for electrons with opposite spin and relative wavevector Q at the same location r.

# **4.2 Inelastic Tunnelling in Spin-Fluctuation Driven Superconductors**

As stated in Sec. 2.5.3, in unconventional superconductors the integrated spin spectrum g 2χ ′′(ω) (see Eq. (2.62)) plays the role of the Eliashberg function. Using the Keldysh diagrammatic formalism for Greens functions, P. Hlobil and J. Schmalian could show that this function is obtainable from the inelastic tunnelling conductance

<sup>3</sup> The wave vector (π/a, π/a) corresponds to h = k = 1/2 in reciprocal lattice units (r.l.u.).

[28]. The inelastic tunnelling conductance due to coupling of electrons to spin fluctuations is given by

$$\sigma^{\rm i}(eV>0) = -\frac{\sigma\_0}{D^2 \nu\_S^0} \int\_{-\infty}^{\infty} \mathrm{d}\omega\_1 \mathrm{d}\omega\_2 \, g^2 \chi''(\omega\_1) \tilde{\nu}\_S(\omega\_2) \times$$

$$\times n\_F'(\omega\_1 + \omega\_2 - eV)[1 - n\_F(\omega\_2)][1 + n\_B(\omega\_1)]\,,\tag{4.1}$$

with the normalized sample DOS ν˜S(ω) = νS(ω)/ν<sup>0</sup> S , Fermi- and Bose-Einstein distribution functions n<sup>F</sup> and nB, as well as the cut-off energy for the bosonic excitation spectrum D. For a two-dimensional Brillouin zone, Eq. (2.62) in the normal state simplifies to

$$
\chi''(\omega) \propto \nu\_S^0 \arctan(\omega/\omega\_{sf}).\tag{4.2}
$$

This leads to σ <sup>i</sup> ∝ g 2 (eV ) <sup>2</sup>/ω<sup>2</sup> sf for small voltages eV ≪ ωsf and σ <sup>i</sup> ∝ g 2 eV /ωsf for eV ≫ ωsf . This behaviour in the normal state is shown in the blue curve of Fig. 4.4(b), which shows the inelastic tunnelling conductance derived from the formal convolution of a constant normal state DOS (blue curve in Fig. 4.4(a)) and an overdamped spin spectrum following Eq. (4.2). In the superconducting state, the spin spectrum is renormalized and a spin gap below ωres evolves (see Sec. 2.5.3). Performing the formal convolution of the superconducting DOS (red curve in Fig. 4.4(a)) and the integrated spin spectrum in the superconducting state (Fig. 2.8(b)), as instructed in Eq. (4.1), one finds the total tunnelling conductances σ tot = σ el + σ i shown in Fig. 4.4(c) for different weights 1/D of the inelastic tunnelling contribution. While the pure elastic tunnelling conductance shows a peak-like feature at ∆ + ωres

**Figure 4.4: Elastic and Inelastic Contributions:** Calculated elastic (a) and inelastic (b) contributions to the differential conductance in superconducting (red) and normal state (blue). The total conductance (c) develops a dip-hump structure around ∆+ωres for high inelastic current contributions (red curves). For small to zero inelastic current contributions (orange/yellow curves) the conductance in the superconducting state stays higher than the conductance in the normal state for eU > ∆. Reprinted figure with permission from [28]. Copyright (2017) by the American Physical Society.

due to coupling to virtual bosons, the inclusion of inelastic tunnelling events leads to a peak-dip-hump structure in the total conductance for energies larger than ∆. The hump at ∼ ∆ + ωres that brings the total conductance in the superconducting state back to match the one in the normal state is due to the sudden onset of inelastic tunnelling events involving real bosonic excitations.

From a series of INS experiments on various cuprate superconductors, the universal relationship between magnetic resonance and superconducting gap of

$$\frac{\omega\_{\text{res}}}{\Delta} = (1.28 \pm 0.08) \tag{4.3}$$

was found [130]. While the resonance mode is directly obtained in INS, more model assumptions, like the applicability of the spin-fermion model, have to be put in, to deduce the spin excitation spectrum from STS and in addition to that, the k-space information is lost in tunnelling experiments. In previous attempts at modelling the Eliashberg function, *a priori* knowledge on its shape had to be put in, which biases the result tremendously due to the ill-posedness of the deconvolution problem [131]. Recently, the application of machine learning algorithms on ARPES data proved to be a powerful concept to reverse-model the spin-spectrum, but this happens at the cost of a number of free parameters which cannot be easily mapped onto physical quantities [132, 133]. Instead, we want to pursue the path of a direct deconvolution of scanning tunnelling data, using *a priori* band structure for the normal state model and the inelastic tunnelling theory described above, as input. That this yields insightful results, has already been demonstrated for monolayer FeSe [27]. The goal here is to extend this deconvolution method to the nodal cuprate superconductors.

### **4.3 Separation of Elastic and Inelastic Tunnelling Events**

The single crystals Bi2Sr2CaCu2O8+<sup>δ</sup> and YBa2Cu3O6+x, that were studied here, were grown by Thomas Wolf at the Institute for Quantum Materials and Technologies (IQMT). The Bi2212 sample is slightly underdoped with a T<sup>c</sup> of 82 K. The Y123 sample is from the batch WAX350−1 and is optimally doped with a T<sup>c</sup> of 92 K.

In search of the bosonic mode at ∆ + ωres in low-temperature scanning tunnelling spectra, it is first important to take a look at the inhomogeneity of the samples to decide which spectra are most promising to study. As reported by Fischer et al., Bi2212 tends to show a large inhomegeneity of its DOS in the superconducting state [120]. This can be confirmed in the direct comparison of the conductance inhomegeneity measured on Bi2212 and Y123, shown in Fig. 4.5. The heat maps of the conductance variation δσ/σ¯(x, y) = PeU<sup>t</sup> −eU<sup>t</sup> (σ(eU, x, y) − σ¯(eU))/σ¯(eU) show that it

**Figure 4.5: Conductance Inhomogeneity**: (a,b) Heat maps showing the variation of the differential conductance over a 50 × 50 nm<sup>2</sup> surface area on Bi2212 (a) and Y123 (b). (c,d) Position averaged bias spectra on the grids shown in (a,b) at T = 0.7 K. The higher inhomegeneity of the Bi2212 surface is reflected in both the conductance variation map and the blurred position averaged spectrum. The characteristic dip and hump are marked by blue and red arrows.

is with ≈ 30 % on Bi2212 about three times higher than on Y123. As a consequence, a position averaged spectrum over a 50 × 50 nm<sup>2</sup> area can preserve detailed gap features better for Y123 (Fig. 4.5(d)) than for Bi2212 (Fig. 4.5(c)). Especially the dip-hump (dip marked by blue, hump marked by red arrows in Fig. 4.6(c,d)) feature is still clearly visible in the position averaged spectrum of Y123 at ϵ ≈ 60 meV but is invisible in Bi2212. Therefore, in the case of Bi2212, an average spectrum at one specific location, at which the dip-hump spectral feature was clearly visible, was chosen for this study. For Y123, the position averaged spectrum was chosen.

Fig. 4.6 shows the measured differential conductance and the numerical derivative of the tunnelling conductance on the two sample surfaces. For Bi2212, the normal state at 84 K was measured with the same tip as the superconducting state at 0.7 K. For Y123, only a spectrum in the superconducting state could be obtained because

**Figure 4.6: Symmetrized Spectra:** First and second derivative of the tunnelling current on Bi2212 (a,b) and Y123 (c,d) in the normal (blue) and superconducting state (black).

drastic tip changes complicated the measurement at elevated temperatures. Due to the large inhomegeneity of the gap of Bi2212 on the surface, the spectrum presented here was recorded at one fixed position whereas the spectrum for Y123 is position averaged over an area of 50 × 50 nm<sup>2</sup> . All spectra were symmetrized, normalized and regularized by Gaussian filtering. The symmetrization accounts for the asymmetric background, seen in the spectra (see Fig. 4.5 (c,d)). Most tunnelling experiments so far, revealed a tilted background conductance that is higher in the negative than in the positive bias range due to the influence of doping on the normal state DOS of Bi2212 and Y123 [120]. Since the tilt here is in the other direction, it is most likely caused by a non-constant tip DOS in the applied bias voltage range. The normalization is chosen such that conductance in the normal and superconducting state are the same for E ≫ ∆. This enables a better comparison with theory, where the normal state is effectively studied at 0 K. In reality, the background conductance is reduced at high temperatures due to increased electron-phonon scattering. The regularization smooths the dI/dU curves for a better display of the numerical derivative and prepares it for the deconvolution procedure in the following section.

The spectrum for superconducting Bi2212 in Fig. 4.6(a) shows a single gap with remanent zero-bias conductance due to the dx2−y<sup>2</sup> gap symmetry and relatively smeared coherence peaks hinting at short quasiparticle lifetimes. This is typical for the underdoped regime and may be caused by its proximity to the insulating phase [120]. Outside the gap, a clear dip of the superconducting spectrum below the normal conducting spectrum, followed by a hump reapproaching it, are visible. The V-shaped conductance in the normal state hints towards strong inelastic contributions to the tunnelling current, as was demonstrated in Fig. 4.4. The hump shows as a peak in the second derivative of the tunnelling current that exceeds the curve of the normal state at ≈ 130 mV in Fig. 4.6(b). The relatively round shape of the superconducting gap in Bi2212 is atypical for a classic d-wave superconductor, in which the naive expectation is a V-shaped conductance minimum. As will be shown, the round shape of the gap can, however, be generated without admixture of an s-wave paring term by respecting the anisotropy of the Fermi surface in the normal state.

The spectrum for superconducting Y123 in Fig. 4.6(c) is qualitatively in excellent agreement with previous STM measurements [134] and shows three low-energy features: (i) a superconducting coherence peak at ≈ 25 meV that is less broadened than in Bi2212, (ii) a high-energy shoulder of the coherence peak and (iii) a lowenergy peak at ≈ 10 meV. A good explanation for the high-energy shoulder and sub-gap peak could be that they arise from the proximity-induced superconductivity in BaO planes and CuO chains [135, 136]. This would certainly fit the fact that these states are missing in the Bi-based compounds and that the sub-gap peak shows a direction-dependent dispersion in ARPES data [129]. The hump lies at ≈ 60 meV. Even without a normal conducting spectrum, the V-shaped background conductance, which is in agreement with the predicted inelastic contribution by magnetic scattering in the spin-fermion model, is obvious in the high-energy regime of the superconducting spectrum.

In order to extract the inelastic part of the tunneling spectrum, it is necessary to determine the elastic contribution by fitting a model function that keeps the complexity as low as possible while still capturing the relevant features of the band structure and pairing strength. To account for the presence of pair-breaking processes, a generalized Dynes function was chosen [137]:

$$\begin{split} \nu\_{s}^{\rm el}(\omega) &= \frac{\sigma\_{0}}{n} \int\_{0}^{\pi/2} \mathrm{d}\varphi \nu\_{0} \nu\_{\mathrm{n}}^{\rm F}(\varphi) \times \\ & \quad \times \left| \mathcal{R} \left( \frac{\omega + i\Gamma(\varphi)}{\sqrt{(\omega + i\Gamma(\varphi))^{2} - \Delta(\varphi)^{2}}} \right) \right| . \end{split} \tag{4.4}$$

Here, n = R dφνn(k<sup>F</sup> , φ) is a normalization factor, ∆(φ) = ∆<sup>0</sup> cos(2φ) is the d-wave pairing potential, Γ(φ) = γ · |∆(φ)| the scattering rate and ν F n (φ) is the angle dependent density of states at the Fermi energy, in the normal conducting phase. The function ν F n (φ) weighs gap distributions for different (kx, ky) by their abundance along the Fermi surface. In order to determine this function, the in plane dispersion is modeled by a tight binding model of the form

$$\begin{aligned} \epsilon(k\_x, k\_y) &= \frac{t\_1}{2} (\cos(k\_x) + \cos(k\_y)) + t\_2 \cos(k\_x) \cos(k\_y) \\ &+ \frac{t\_3}{2} (\cos(2k\_x) + \cos(2k\_y)) + \frac{t\_4}{2} (\cos(2k\_x) \cos(k\_y) \\ &+ \cos(k\_x) \cos(2k\_y)) + t\_5 \cos(2k\_x) \cos(2k\_y) \\ &- \mu, \end{aligned} \tag{4.5}$$

with chemical potential µ and hopping parameters t as proposed in Ref. [138] for a near optimally doped Bi2212 crystal. For Y123, the model function [139]

$$\epsilon(k\_x, k\_y) = 2t\_1 \cos(k\_y) 2t\_3 (\cos(2k\_x) + \cos(2k\_y))\tag{4.6}$$

$$\pm \left[ 4 \cos^2(k\_x) (t\_1 + 2t\_2 \cos(k\_y))^2 + V^2 / 4 \right]^{1/2} + \mu \tag{4.7}$$

was used, where the fourth term accounts for the ortho-II band folding due to CuO chains. The parameters t<sup>1</sup> = −0.558 eV, t<sup>2</sup> = −0.31t1, t<sup>3</sup> = −0.5t2, µ = 0.17 eV and V = 0.025 eV were used. The parameters µ and V were adjusted such that the three resonances in the DOS from the quasi-1D band correspond to the energetic positions of peaks in the spectrum. Compared to Ref. [139], the lower µ corresponds to lower hole-doping and shifts the lowest energy resonance to lower energies compared to the high energy resonance (coherence peak shoulder). A reduced ortho-II potential V leads to a more pronounced low-energy resonance due to the LDOS peak in the near-nodal direction as seen in Fig. 4.7 (d).

An analytic expression for the Fermi wave vector k<sup>F</sup> (φ) is retrieved from the solution of ϵ(k, φ) − µ = 0 where ϵ(k, φ) is the polar representation of Eq. (4.5). The normal DOS is then given by

$$\nu\_{\mathbf{n}}^{\mathbf{F}} = \oint \frac{\mathbf{d}^2 k}{|\nabla\_{\mathbf{k}} \epsilon|} \to \int\_l \frac{\mathbf{d}\varphi}{\nabla\_{k,\varphi} \epsilon(l(\varphi))} \left\| l'(\varphi) \right\|\_2 \,, \quad \varphi \in M\_{\varphi}, \tag{4.8}$$

where l(φ) = (k<sup>F</sup> (φ) cos(φ), k<sup>F</sup> (φ) sin(φ))<sup>T</sup> is a parametrization of the path along the Fermi surface, ∥·∥<sup>2</sup> is the Euclidean norm and M<sup>φ</sup> = {φ|k<sup>F</sup> (φ) ∈ 1.BZ}.

The calculated Fermi surface contour and DOS in the normal state are shown in Fig. 4.7. The fits to the superconducting spectra using Eq. (4.4) and parameters listed in Tab. 4.1 are shown in Fig 4.8(a,b). Unlike in fully gapped superconductors, the angle-integrated inelastic spectrum can be non-zero down to vanishing bias voltage because ∆(k) has a nodal structure. This prevents a direct assignment of the differential conductivity for e|U| ≲ ∆ to the purely elastic tunneling contribution like it is possible for the s<sup>±</sup> superconductor monolayer FeSe [27, 28]. Instead, the fact that

**Figure 4.7: Normal State Electrons:** (a,c) Calculated Fermi surface in the first 2D BZ of Bi2212 (a) and Y123 (c). For Y123 the contribution of the quasi-1D ortho-II band is shown in a dashed blue line. (b,d) Calculated Fermi wavevector k<sup>F</sup> (φ) (orange) and normal conducting DOS along the Fermi surface contour ν F <sup>n</sup> (φ) (purple) as function of polar angle in the first BZ quadrant for Bi2212 (b) and Y123 (d). For the ortho-II band, the Fermi wavevector and normal DOS are drawn as blue and green dashed lines in (d).

σ tot = σ el + σ inel and the physical constraints σ inel(0) = 0 and σ inel(e|U| > 0) > 0, shall be used. As demonstrated in Sec. 4.2, σ inel is in principle a convolution of ν<sup>s</sup> and the bosonic spectral function. Hence, for a slowly varying bosonic function, which is expected in the range 0 < e|U| < ∆, σ inel is dominated by the elastic contribution νs(ω) with some factor, say (1 − α). Consequently, in the range 0 < e|U| < ∆, σ el is well guessed by our Dynes fit times a factor α < 1. For Bi2212 α is chosen in the following way: From the requirement that the integrated elastic contribution to the conductance in the spectroscopic range is equal for normal and superconducting state, one can find the scaling factor for the elastic part that best describes its proportion to inelastic processes. This is shown in Fig. 4.9(a) in Sec. 4.4. For Bi2212 this proportionality factor is α = 0.6. In order to keep the condition σ el(∞) = σ0, the experimental curve is scaled up by 1/α instead of scaling down the fitted curve. The resulting elastic and inelastic contributions are shown in red/grey in Fig. 4.8(c,d). For Y123, only the condition σ inel > 0 was enforced by using the scaling factor α = 0.9 because normal state spectra were lacking. The influence of this scaling factor on the extracted bosonic spectrum is discussed in detail in the next section.

**Table 4.1: Fit Parameters**: Fitted Dynes function (Eq. (4.4)) parameters for Y123 and Bi2212 superconducting spectra.


### **4.4 Determining the Dispersion of the Bosonic Glue**

**Figure 4.8: Spin Spectrum Extraction:** (a,b) Elastic conductance fit for eV ≲ ∆(red) following Eq. (4.4) with parameters from Tab. 4.1 to the measured conductance on Bi2212 (a) and Y123 (b). (c,d) Top: The measured conductance spectrum is scaled up by 1/α. The inelastic part is the difference between total conductance (black line) and the elastic part (red line). Red and orange dashed lines show the total conductance obtained from convolving the elastic conductance guess with the orange and green spin spectra. Bottom: Spin spectra obtained from direct deconvolution in Fourier space (orange) and from the Gold deconvolution algorithm (green). Results for Bi2212 are shown in (c), for Y123 in (d).

In order to extract the spin spectrum, Eq. (4.1) at T = 0 is rewritten as a formal convolution integral:

$$\sigma^{\rm inel}(eU) \propto \int\_0^{eU} \mathrm{d}\omega \nu\_s(\omega) g^2 \chi''(eU - \omega)$$

$$= \left( [\nu\_s \cdot \Theta] \* \left[ g^2 \chi'' \cdot \Theta \right] \right)(eU) \,. \tag{4.9}$$

As can be seen from Eq. (4.9), the integrated spin spectrum χ ′′ is in information processing terms the source, the DOS in the superconducting state is the kernel and the inelastic tunneling conductivity is the signal. Θ denotes the Heaviside step function. In general, retrieving the function g 2χ ′′ from the inelastic conductivity is an ill-posed problem. Additional aspects that complicate the problem are the following:


We compare two methods with which the Eliashberg function was determined: By direct deconvolution in Fourier space and by the Gold algorithm [140, 141]. The advantage of the Gold algorithm is that for a positive kernel and signal, the result of this iterative deconvolution method is always positive, in agreement with our physical constraint. It is explained in more detail in Appendix A.

The regularized bosonic function from direct deconvolution in Fourier space

$$g^2 \chi''(\omega > 0) = \mathcal{F}^{-1}\left(\frac{\mathcal{F}(\sigma^{\text{inel}})(t)}{A\sqrt{2\pi}\mathcal{F}(\nu\_s \cdot \Theta)(t)}\right) \tag{4.10}$$

is shown in Fig. 4.8(c,d) in orange. The contributions at low energies mostly stem from the scaling factor α. The negative part of the spectrum for E < ∆ is due to the abrupt change in elastic conductivity at zero energy (multiplication with Heaviside distribution) which results in heavy oscillations in time space. Therefore, the circular deconvolution contains nonzero contributions for E < 0 and negative contributions for 0 < E < ∆ which are not physical. Despite these shortcomings, the bosonic spectrum recovers well the tendency of the total conductivity (Fig. 4.8(c,d) dashed orange) and shows the expected behaviour at medium and high energies, i.e. a resonance at ∆ < E < 2∆ and approach of the normal state bosonic function for E < 3∆.

Using the result of the direct deconvolution method as a first guess to the Gold algorithm, the bosonic spectrum shown in green in Fig. 4.8(c,d) was obtained after 2,000,000 iterations. Again, the high value at E = 0 is a consequence of the scaling with factor α. Negative contributions are now gone and the bosonic function recovers well the total elastic conductivity shown in Fig. 4.8(c) (dashed green), especially the dip-hump structure, and shows a very clear resonance at ωres ≈ 70 meV ≈ 1.1∆<sup>0</sup> ≈ 1.2∆max ≈ 1.5∆¯ for Bi2212 and at ωres ≈ 60 meV ≈ 2.0∆<sup>0</sup> ≈ 2.4∆max ≈ 3.6∆¯ for Y123.

The resonance mode of Bi2212 extracted in this work is higher in energy than reported in inelastic neutron scattering experiments (ωres ≈ 43 meV at the AFV) [26] and closer to the resonance determined by optical scattering (ωres ≈ 60 meV) [142]. Due to the loss of k-space information in tunnelling, the centre of the resonance is expected to be shifted to higher energies compared to the INS results (compare Fig. 2.8(a) and (b)). The ratio ωres/∆max lies within the current range of error of Eq. (4.3). In most other extraction methods of the bosonic mode energy, the normal state DOS is not respected, which is why the ∆<sup>0</sup> used there is most similar to what is here called ∆max. ∆max is the largest gap value that contributes to the elastic conductance spectrum.

For Y123, the resonance mode is significantly higher in energy in the bosonic density spectrum than experimentally found by INS in optimally doped samples with ωres ≈ 41 − 46 meV [121, 130] and even lies outside the onset of the spin scattering continuum at ω<sup>c</sup> ≈ 59 meV in the slightly underdoped regime [130]. Apart from the k-space integration, which shifts the peak centre to higher energies, several other factors can play a key role: (i) The well-studied 41 meV odd-parity mode is paired with an even-parity mode at ω e res ≈ 53 − 55 meV [130, 143–145] which may be of the same origin as it vanishes at Tc. This mode appears with a ≈ 3 − 20 times lower intensity in INS than the odd-parity mode, but this does not necessarily have to hold for a tunnelling experiment. (ii) The bosonic spectrum extracted here is essentially poisoned by phononic contributions from every k-space angle. Even INS data at specific q-vectors have to be corrected by known phonon modes. In the integrated spectrum obtained here, this is simply not feasible anymore. (iii) Apart from physical arguments, there can also be made sceptical remarks on the

**Figure 4.9: Numerical Scaling Factor** α**:** (a) Determination of α happens through the boundary condition R dϵσel <sup>n</sup> = R dϵσel s . (b) Variation of α leaves the general shape of the extracted bosonic spectrum unaffected except for the magnitude of its zero-energy peak.

deconvolution procedure: Evidently, it heavily depends on the guess of the elastic tunnelling conductance, which in this case does not contain strong-coupling features from an Eliashberg theory, like in Fig. 4.4. The large deviation from ωres/∆ ≈ 1.28 is most likely caused by the determination of ∆. For gap sizes, the authors from Ref. [130] use a constant ∆max = 39.5 meV throughout a large range of doping in the underdoped regime. In a collection of STS results spanning over 20 experiments, however, the extracted gap value is lying between ∆ = 20 − 30 meV for underdoped and optimally doped samples [120]. Even with inclusion of the copper oxide chain gap the spectroscopic results on Y123 in this work do not support an effective gap value of > 30 meV because the total conductivity is already on the decrease at this energy.

In order to make sure that the introduction of the numerical scaling factor α has no poisoning effect on our extracted spin spectrum, the deconvolution of the Bi2212 spectrum by Gold's algorithm was performed for four different values of α. The results shown in Fig. 4.9 are comforting in the sense that the overall shape of the spin spectrum is unchanged. The only major difference lies in the magnitude of the zero-energy peak which is to be expected from a scalar multiplication, but since this peak is anyhow out of the bounds of physical contributions it does not harm the analysis. The oscillations in the obtained function g 2χ ′′ for high energies are also not affected by a change of α. They are caused by the FFT of the elastic conductance, so their position and spacing only changes upon changing ∆. They are a side effect of the algorithm that might be reduced by clever filtering of the Fourier transformed elastic conductance that enters the algorithm.

### **4.5 Summary**

The STM junction geometry was utilised to effectively couple to real bosons in inelastic tunnelling spectroscopy (IETS) on cuprate superconductors. The superconducting and normal state tunnelling spectra showed significant inelastic contributions for underdopedBi2212 and optimally dopedY123. Since the effective Eliashberg function is inaccessible through the normal state in unconventional superconductors, it was obtained by deconvolution of the superconducting spectrum. In contrast to fully gapped unconventional superconductors like monolayer FeSe [27], the separation of elastic and inelastic tunnelling contributions proves to be more difficult in the nodal superconductors. Particularly, the k-dependence of ∆ and the normal state DOS are important factors that were respected in a tight-binding model of the superconducting DOS based on ARPES measurements. By using the Gold algorithm for the deconvolution procedure, a strictly positive effective Eliashberg function was obtained. It takes the form of the expected integrated spin spectrum in the spin-fermion model with a resonance mode at ωres > ∆0. The resonance mode was found at ωres ≈ 70 meV in Bi2212 and at ωres ≈ 60 meV in Y123. Despite the loss of momentum information in STS, a description in momentum averaged electronic and bosonic spectral functions yields reasonable results because the tunnelling spectrum is dominated by the hot spots on the Fermi surface.

# **5 Yu-Shiba-Rusinov States and Long-livedBosonic Excitations in Granular Aluminium**

Granular aluminium (grAl) has recently establisheditselfin parts of the superconducting qubit community as a suitable superinductor in Al/AlO<sup>x</sup> based superconducting circuits [29–31] and as a template for magnetic-field resilient qubit fabrication [32, 33]. Not only are excess quasiparticles identified as a limiting loss mechanism in these circuits, transport measurements in grAl near its MIT also hint at localized magnetic moments that might be responsible for the 1/f flux noise in Al/AlO<sup>x</sup> based junctions [34, 35]. STM provides a local probe of the quasiparticles on a microscopic scale and is sensitive to spin scattering, which should manifest as a pair of Yu-Shiba-Rusinov (YSR) states below Tc. In this chapter, the first low-temperature STM/STS results on grAl samples in two different resistivity regimes are presented and put into perspective by comparison with a polycrystalline pure aluminium film. The main focus in this work lies on the sub-gap YSR states that were found in the high resistivity grAl samples together with the inelastic excitations that are likely connected to phase modes of the superconducting order parameter. Lastly, an experimental procedure to systematically locate the spots, at which localized magnetic moments reside, is outlined and first results are shown.

**Disclaimer**: I, the author, want to emphasize that I was not personally involved in the measurements of the grAl samples but performed the measurements of the polycrystalline Al films (Chapter 5.2.1, 5.7). The measurements of the grAl samples were performed by the former group member Dr. Fang Yang and his at the time master student Tim Storbeck in the group of Prof. Wulfhekel at Karlsruhe Institute of Technology (KIT). I was involved in data analysis and presentation, as well as manuscript preparation for the publication that followed in Physical Review B. Parts of the results presented in this thesis chapter have been published by me as a co-author in [146]. Reused data is marked by the self-citation in the respective figure caption.

## **5.1 Relation between** T<sup>c</sup> **and Oxygen Partial Pressure**

The discovery of granular superconductors and their enhanced superconducting transition temperature, compared to the homogeneous bulk metal, has its roots in the group of Ben Abeles in Princeton [20, 147]. In the following years, grAl quickly became the centre of attention due to its controlled growth [148] and the large T<sup>c</sup> enhancement to up to 3.15 K [149]. The term "granular superconductor" refers to a system in which superconducting grains are separated by insulating barriers, effectively building an array of intercoupled Josephson junctions (JJ). For aluminium, this is realized by evaporating aluminium onto a room temperature (or colder) substrate in an oxygen atmosphere. At oxygen partial pressures larger than roughly 1×10<sup>−</sup><sup>5</sup> mbar [148], the aluminium partially oxidizes in such a way that pure aluminium grains of a well-defined grain size a are separated by a non-stoichiometric AlO<sup>x</sup> matrix. The normal state DC resistivity ρdc is directly proportional to the oxygen partial pressure and is the parameter of choice to compare samples grown in different setups. Up to the Mott resistivity ρ<sup>M</sup> ≈ 400 µΩ cm, the mean size of the Al grains on room temperature substrates gradually shrinks to a = (3 ± 1) nm and further remains constant for larger resistivities, yet with a smaller variance [150]. Close to the Mott resistivity, T<sup>c</sup> peaks [151, 152] and reaches 2.17 K (for room temperature substrates) [153] before quickly decreasing for samples with ρdc ≳ 3 × 10<sup>3</sup> µΩ cm due to a superconductor-to-insulator transition (SIT) [149]. This relation between normal state resistivity and T<sup>c</sup> creates a superconducting dome, which is depicted in the phase diagram in Fig. 5.1. According to Ref. [153], the superconducting gap size ∆ follows the general trend of T<sup>c</sup> in its dependence on DC resistivity below ρ<sup>M</sup> , but jumps near the Mott resistivity, which leads to a change of the BCS coupling strength ∆/T<sup>c</sup> from ∼ 1.78 to ∼ 2.1 for the high resistivity regime [153]. The maximum gap size determined in their THz absorption experiments was ∆ = (0.38 ± 0.04) meV for a sample with ρdc = 1600 µΩ cm.

#### **Figure 5.1: Superconducting Dome**:

The superconducting transition temperature Tc of grAl (black line) gives rise to a superconducting dome in the phase diagram when plotted over the logarithmic normal state DC resistivity ρdc. It peaks at the Mott resistivity ρ<sup>M</sup> ≈ 400 µΩ cm, at which the material features a MIT for T > Tc. Sketched according to data from Ref. [153].

Until today, no general consensus has been reached on the underlying mechanisms that govern the T<sup>c</sup> enhancement, drive the SIT and allow superconductivity to survive for resistivities far above the Mott limit. A few widely accepted concepts shall, however, be briefly mentioned here. For further reading, the reader is referred to the correspondingly cited literature. Suggested mechanisms include

	- *Size effects* of the discrete energy spectrum inside the grain lead to an effective BCS-like electron-electron interaction inversely proportional to the grain size and a transition to a strong coupling limit for small grain sizes [154].
	- *Shell effects* in the delocalized states of metallic nanosized clusters containing a specific number of atoms can lead to an increased pairing strength [155, 156].
	- *Phonon softening* due to the reduced dimension of the superconductor leads to a down shift of phonon DOS in energy and an increased Eliashberg electron-phonon coupling [157, 158].

All these effects focus on an enhanced superconductivity of the metal grains. In the presence of the insulating barriers (AlOx) they have to be combined with an intergrain coupling mechanism, e.g. a Josephson coupling like proposed in Ref. [154, 159, 160], in order to explain an increase in bulk Tc.

	- *Quenching of superconductivity in individual grains* happens due to large amplitude fluctuations of the order parameter [161, 162] when the Josephson coupling between them is weak.
	- *Loss of macroscopic phase coherence* is caused by the systems increasing susceptibility to phase fluctuations of the order parameter (Experimental studies) [149, 153].
	- *Percolation* arguments can be made for the superconducting junctions not breaking down all at once but in a statistical manner leaving room for a few superconducting channels via well-coupled grains [160] (This can easily be discarded if a local STM probe of the LDOS shows a superconducting gap on each grain).

# **5.2 Microscopic Grain Coupling and the Superconducting Gap**

The aluminium and granular aluminium films were grown on Nbdoped (0.7 weight %) SrTiO<sup>3</sup> (Nb:STO) single crystals of (001) orientation. The substrate provides a conducting, flat and nonreactive surface. Aluminium was deposited by electron beam evaporation in a commercial PLASSYS MEB550 at a rate of 1 nms<sup>−</sup><sup>1</sup> up to the desired film thickness and the grAl films were grown in an oxygen atmosphere. The growth parameters and estimated sample resistivities are summarized in Tab. 5.1. Subsequently, the samples were immediately transferred to the STM setup via a vacuum suitcase to keep surface contamination to a minimum.

**Figure 5.2: Sample Structure**: Sample constituents and schematic structures: Left: The grAl samples feature nano-sized aluminium grains embedded in an amorphous AlOx matrix. The coupling of individual superconducting Al grains essentially leads to a large array of Josephson junctions. Right: The polycrystalline Al film features larger crystallites and macroscopic phase coherence of the superconducting order parameter.

**Table 5.1: (Granular) Aluminium Samples**: Investigated samples including their estimated resistivity ρdc at oxygen partial pressure pO<sup>2</sup> and film thickness d following established growth recipes further detailed in Ref. [29].


Fig. 5.2 schematically summarizes the most important structural properties that are expected from the grAl and Al films. Based on early electron microscopy studies [148, 150], the grAl film should consist of nano-sized Al grains separated by an amorphous AlO<sup>x</sup> matrix. The thickness of the oxide, which is fine-tuned by the oxygen partial pressure during growth, eventually causes the SIT. As long as Cooper pairs can tunnel between two neighbouring Al grains i and j with order parameter ∆eiφi,j , the material is superconducting and acts as a large array of Josephson junctions with individual Josephson couplings Jij . Above a critical thickness the grains uncouple and the material becomes insulating. The polycrystalline Al film should consist of larger crystallites due to the high diffusion mobility of Al at room temperature [146]. There, the phase of the superconducting order parameter is expected to be universal throughout the sample.

#### **5.2.1 Pure Al Film**

The measurement of sample A, the pure Al film, served as a control experiment. Fig. 5.3(a) shows the topographic image of the film that features crystallites of about 50 nm lateral size. The large-range I(U) curve in Fig. 5.3(b) confirms the metallic behaviour for large bias voltages. Fig. 5.3(c) shows the low bias dI/dU spectrum at T = 30 mK which reveals a superconducting gap of ∆ = (191 ± 1) µeV. The gap size was determined by a fit to the BCS DOS. This is close to the value of 200 µeV reported for a 42 nm thick Al film [170]. In all structural and electronic aspects, the sample behaves as expected.

### **5.2.2 Oxygen Poor GrAl Film**

Next, the oxygen poorgrAl sample G1 with an estimated resistivity of ρdc = 300 µΩ cm was investigated. As can be seen in the topographic image in Fig. 5.4, a granular structure of round Al grains with sizes ranging from 5 − 10 nm in diameter formed.

**Figure 5.3: Pure Al Sample**: (a) Topographic image of the pure Al film (U = −1 V, I = 100 pA) showing crystallites of about 50 nm lateral size. (b) Typical large range I(U) curve displaying ohmic behaviour of the sample. (c) dI/dU spectrum at T = 30 mK showing the superconducting gap of ∆ = (191 ± 1) µeV. The gap size was determined from the fit of a BCS gap (blue line). Data previously published in [146].

For objects this small and pronounced, the shape of the tip plays an important role in the topographic image [171]. In fact, the resultant topography in the constant current mode reflects a convolution of tip surface and sample surface. The tip cannot be etched to arbitrary sharpness and tip apex radii of a few nm are not unusual [112, 172]. Taking this into account, the observation agrees well with previously reported grain sizes and spread for similarly grown grAl films [150]. A comparison between typical I(U) curves measured on different grains of G1 (Fig. 5.4(b)) and the I(U) characteristic of the pure Al film A (Fig. 5.3(b)) confirms that sample G1 is still in the metallic regime.

Fig. 5.4(c) shows the dI/dU spectrum at T = 30 mK with a superconducting gap of ∆ = (298 ± 1) µeV, which is significantly larger than for the pure Al film. The gap size was determined by a BCS fit (blue line). The gap size was practically identical when probed on different grains (see Fig. 5.4(d)). For similar grAl films, the authors of Ref. [173] derived ∆ = 344 µeV from Tc. THz spectroscopy experiments yielded ∆ = (336 ± 8) µeV [153]. The uniform gap size on the surface for different grains could be explained by either a homogeneous gap enhancement mechanism which is independent of the grain size or a strongly coupled network in which the ∆ of different grains are aligned. Either way, grains on the surface would exhibit

**Figure 5.4: Oxygen Poor Sample**: (a) Topographic image of the grAl film (U = 60 mV, I = 240 pA) with ρdc ≈ 300 µΩ cm showing grains of 5−10 nm lateral size. (b) Typical large range I(U) curves displaying metallic behaviour of the sample. (c) dI/dU spectrum at T = 30 mK showing the superconducting gap of ∆ = (298 ± 1) µeV. The gap size was determined from the fit of a BCS gap (blue line). (d) dI/dU spectra on different grains reveal a practically universal gap except for the grain marked in cyan in (a), which exhibits an additional pseudogap of roughly 8 ∆. Data previously published in [146].

the same gap as in the bulk making the STS measurements here comparable to bulk measurements. In Ref. [173], where ∆ was derived from T<sup>c</sup> measurements, an overestimation can be explained by the overall higher T<sup>c</sup> compared to Ref. [153], which indicates a significant difference in the systematic approach. It should be noted that THz spectroscopy represents a real bulk measurement of ∆, whereas in STS the sample LDOS at the tip position is measured. Other factors, that come into play, are: 1. the overall lower oxygen concentration near the surface for the UHV samples used in this STM study compared to the ambient samples from different studies, 2. the overall higher oxygen concentration at the surface compared to the bulk due to the growth procedure, 3. the fact that only a small area on the surface was probed by STM and 4. the limited comparability of samples produced by different recipes and in slightly different environments. Considering all these factors, the agreement between the gap size determined here and bulk measurements is surprisingly good.

**Figure 5.5: Josephson Tunnelling**: dI/dU spectrum under extreme tunnelling condition (high tunnelling rate) on the insulating grain (a) from Fig. 5.4(a) and a more metallic grain (b). Tip stabilization condition: U = 500 µV, I = 200 pA. Data previously published in [146].

On the grain marked in cyan in Fig. 5.4(a), the low energy spectrum shows an additional pseudogap<sup>1</sup> of roughly 8 ∆, as shown as by the cyan line in Fig. 5.4(d). As we will see, this behaviour is rather atypical for the oxygen poor sample but not uncommon for the oxygen rich sample. The LDOS for eU > ∆ is heavily reduced in comparison to the other grains in this sample. For extreme tunnelling conditions on this insulating grain, i.e. a conductance at the scale of µS and a tunnelling current approaching the critical current, the dI/dU spectrum shown in Fig. 5.5 (a) was obtained. It features a sharp peak at zero bias, flanked symmetrically by in-gap peaks and a dip to negative differential conductance. The spectrum for a metallic grain at the same tunnelling conditions does not show these features (Fig. 5.5(b)). Under extreme tunnelling conditions, the voltage drop over the vacuum barrier seems to become comparable to the voltage drop over the sample. As a consequence, the measured dI/dU spectrum also contains information on the electronic transport inside the sample. The zero-bias peak is reminiscent of the Josephson DC current in an SIS junction, which can be understood as follows: A Cooper pair in one superconductor is destroyed through an Andreev reflection at the insulator interface and continues as an evanescent wave in the insulator. If the evanescent wave overlaps with the one coming from the opposite superconductor, then a Cooper pair can be created in the opposite superconductor by an additional Andreev reflection of the electron-hole pair. Hence, a dissipationless DC current flows through the superconductor. The oscillations and negative differential conductance at e|U| ≳ 0 are consistent with previous STS measurements of the Josephson DC current and explained within the P(E) theory [106, 175]. Direct lead-to-lead elastic cotunnelling and non-local Andreev reflections [176, 177], which cause very similar dI/dU features [178], are

<sup>1</sup> Although the term "pseudogap" is nowadays often associated with high-temperature superconductivity, it is here meant in the original sense popularized by Mott, i.e. a depletion of electronic states near the Fermi energy due to increased localization of electronic wave functions near an MIT [174].

not expected because the superconducting coherence length in grAl is reduced to about 10 nm [179], so below the film thickness. If the additional features near zero bias voltage stem from intergrain tunnelling processes, then the dI/dU spectrum should also exhibit coherence peaks at 2∆, due to the doubled energy needed to transfer a quasiparticle from one superconductor to another superconductor with the energy gap ∆. Unfortunately, no data in this bias range was recorded. It leaves open the question, whether the acquired dI/dU spectrum is the superposition of the LDOS on the grain and intergrain tunnelling characteristics in the bulk or a different phenomenon.

#### **5.2.3 Oxygen Rich GrAl Film**

In Fig. 5.6, the general behaviour of the oxygen rich film with ρdc = 2000 µΩ cm is summarized. The topographic image in Fig. 5.6(a) confirms that its grain sizes are similar to the oxygen poor film's, an observation well in agreement with previous electron microscopy studies [150]. The much higher resistance of this sample is, however, reflected in the fact that I(U) curves (Fig. 5.6(b)) measured on individual grains, marked in Fig. 5.6(a), exhibit insulating behaviour, so a region of almost zero current for voltages lower than a threshold of roughly 70 to 300 mV (depending on the grain and bias sign). Above the threshold voltage, a differential conductance (I(U) slope) comparable to the oxygen poor sample is obtained. The sudden onset of electronic transport above the threshold voltage signals that the grains are in the Coulomb blockade regime. At smaller tunnelling rates (lower tunnelling conductance) the I(U) spectrum evolves to a Coulomb staircase, which is reflected in the dI/dU spectrum in 5.6(c) as resonant tunnelling peaks for e|U| > Ec, where E<sup>c</sup> is the charging energy [180].

**Figure 5.6: Oxygen Rich Sample**: (a) Topographic image of the grAl film (U = 400 mV, I = 180 pA) with ρdc ≈ 2 × 10<sup>3</sup> µΩ cm showing grains of 2 − 5 nm lateral size. (b) Typical large range I(U) curves obtained at the grains indicated by the same colors in (a). The I(U) curves display insulating behaviour of the sample with varying gaps. (c) dI/dU spectrum of a particularly isolated grain in the Coulomb blockade regime showing charging peaks at the energetic positions indicated by triangular markers. Data previously published in [146].

**Figure 5.7: Grain Coupling**: (a-f) Constant height maps of the area shown in Fig. 5.6(a) at the indicated bias voltages U display the local tunnelling current I. Multiple neighbouring grains form electrical clusters that charge and discharge collectively. Data previously published in [146].

The energy of these peaks corresponds to first, second, etc. charging state of the single grain or cluster [181]. Using

$$E\_c = \frac{e^2}{2C},\tag{5.1}$$

the first charging state at U ≈ ±0.2 V in Fig. 5.6(c) corresponds to a capacitance of C ≈ 0.4 aF. The threshold voltages or charging energies ranging from ±70 to ±300 meV thus correspond to capacities between 0.3 and 1.1 aF. Approximating the capacitor as a sphere, of which one half is embedded in a dielectric of permittivity ε<sup>r</sup> = 9 and thickness 0.5 nm, the smallest capacitance corresponds to a charged object of approximately 2 nm in diameter (see Sec. 5.5 for the model). Apparently, the smallest chargeable object is a single grain, but larger clusters of grains can also form electrical islands that are simultaneously charged and discharged. This conclusion is supported by the tunnelling current maps shown in Fig. 5.7. The maps were recorded in the constant-height-mode on the area from Fig. 5.6(a). An abrupt increase in tunnelling current for increasing bias voltage is seen at once for several neighbouring grains, indicating that several grains indeed form strongly coupled electrical clusters. Recent measurements of a transmon qubit made from thin film grAl also find evidence for a strong electrical coupling of about 10 grains to one effective Josephson junction [32].

The detailed discussion of the low energy spectrum and especially the superconducting gap in sample G2 is shifted to the sections 5.3 and 5.4, which focus on the in-gap and out-gap structure of the DOS respectively. It should only be aforementioned here, that the gap size for this sample is also enhanced by the granularity and with ∆ = (312 ± 1) µeV also slightly larger than for sample G1, which is in agreement with optical spectroscopy measurements [153] on samples of similar resistivity.

### **5.3 Spectroscopic Signs of YSR states**

Fig. 5.8(a) shows a typical high-resolution spectrum of the superconducting gap in sample G2. A BCS fit (blue line) yields a superconducting gap of ∆ = (312 ± 1) µeV. Remarkable in this spectrum are the two in-gap states at energetically symmetric positions with respect to the Fermi energy, which are indicated by black triangles. These in-gap states are found for eleven out of 23 individually probed grains. The intensity of these peaks is low compared to the differential conductance outside the gap. Their energetic position is very close to ∆ and thus rules out that they are resonances due to Andreev bound states (cf. Fig. 5.5). The observation of the Kondo effect in grAl films for temperatures larger than T<sup>c</sup> [35, 163] indicates that the resonances seen here are spectroscopic signs of Yu-Shiba-Rusinov (YSR) states. In simplest terms they can be thought of as the analogon to the Kondo effect for a superconductor. They are caused by the exchange coupling between conduction electrons and unpaired localized electrons or magnetic moments in the presence of

**Figure 5.8: YSR States**: (a) dI/dU spectrum on an individual grain of the oxygen rich sample. Triangular markers indicate the positions of YSR states inside the gap. (b) Heatmap of the dI/dU spectrum as a function of lateral tip position across a single grain. The line profile is indicated in orange on top of the grey-scale topographic image. YSR states move symmetrically in energy as a function of the tip-grain distance. Data previously published in [146].

a superconducting gap. In this scenario, a Cooper pair can be broken for energies smaller than ∆ if the resulting pair of quasiparticles participates in the screening of the localized spins and thus compensates the cost of kinetic energy by a reduction of the mean potential energy.

As pointed out by the authors of Ref. [35], two scenarios for the origin of unpaired magnetic moments have to be distinguished: They can either reside on the Al grains or in the insulating oxide. In the former scenario, a sufficient decoupling of individual grains in the Coulomb blockade regime can result in an odd number of electrons on individual grains in the equilibrium state. These so-called Kubo spins are suggested as the localized moments responsible for the Kondo effect at T > T<sup>c</sup> by the authors of Ref. [182]. It still remains unclear, how exactly these localized and unpaired electrons would survive the onset of superconductivity, where electrons in the nanosized Al grains pair up to Cooper pairs. The quenching of superconductivity in a few grains by a combination of size effects and decoupling from the rest of the network at first glance contradicts the observed large T<sup>c</sup> and enhanced electron-phonon coupling [153] for ρ > ρ<sup>M</sup> . In the latter scenario, the nonstoichiometric composition of AlO<sup>x</sup> could lead to trapped, unpaired electrons in dangling bonds at surfaces or interfaces. The existence of these magnetic impurities would not be altered by the onset of superconductivity and their abundance would scale with the thickness of the oxide, so the resistivity of the samples.

In order to figure out, which scenario is present, the dependence of the YSR states on an electric field was measured. This was realized by moving the tip laterally across a single grain and recording the full dI/dU spectrum in equidistant steps. The line profile and spectroscopic positions with respect to the grain are shown in Fig. 5.8(b) in the top panel. In the lower panel the dI/dU spectrum is shown in the form of a heatmap. While the superconducting gap edge (dark red) remains fixed, the YSR states (apricot colour) vary in energy as a function of tip position. If the unpaired spins resided in electrically decoupled grains with an odd number of electrons, then the electric field of the biased tip-sample junction would shift the chemical potential of the grain. As a result, the energetic position of the gap edge would move but the position of the YSR states with respect to ∆ would stay fixed. This contradicts the results shown here. Instead, unpaired spins in the oxide could explain the experimental results. The voltage drop over the insulating oxide barrier between the grains is affected by the electric field at the surface. With a variation of the voltage drop across the insulator, the energy of the unpaired spins is changed and ultimately their exchange coupling to the conduction electrons is modified. As demonstrated by the authors of Ref. [183], this effect leads to an energy shift of the YSR states with respect to ∆ as a continuous function of the exchange coupling.

## **5.4 Spectroscopic Signs of Long-lived Plasma or Higgs modes**

Turning the attention to the out-gap structure of the differential conductance displayed in Fig. 5.9 (black line), which was recorded on sample G2, a pseudogap like in Fig. 5.4(d) of > 8∆ is retrieved<sup>2</sup> . Only now, a previously absent or invisible series of peaks inside this pseudogap is apparent. The energetic positions and width of these peaks are too small to be caused by charging effects, like in Fig. 5.6(c). The sharpness of the peaks is not in agreement with simple pair-breaking by inelastic tunnelling, which would make these excitations strongly damped because the electron pairs quickly recombine. Instead, their shape resembles the coherence peaks. Out-gap copies of the coherence peaks are indicators of inelastic tunnelling processes involving a rather long-lived bosonic mode: The incident tunnelling electron excites a boson of defined energy and the final scattered electron with E ≳ ∆ ends up in the quasiparticle resonance of the BCS single particle DOS. Due to the high DOS of the final state electron, this process has a large weight and prominent copies of the coherence peaks, shifted from their original position by the boson energy, are the consequence. The first copy corresponds to the creation of one boson with energy ω, the second to two bosons, etc.

**Figure 5.9: Inelastic Excitations**: High resolution dI/dU spectrum (black line) featuring a superconducting gap and several inelastic excitations which appear as copies of the coherence peak for eU > ∆. The excitation ladder can be modeled by five energetically symmetric inelastic excitations and fit to the spectrum (red line). The modeled spectrum for zero, one, two, etc. excitations is individually displayed in blue, green, violet, and so on. Except for the excitation shown in violet, all excitation energies are integer multiples of ωexc. The energies ∆ + n · ωexc are indicated by dark grey dashed lines. The inset reproduces THz spectroscopy measurements of samples with ρ ≈ 2000 µΩ cm from Ref. [153]. Data previously published in [146].

<sup>2</sup> The discussion of this pseudogap is postponed to the end of this section.

**Table 5.2: Fit Parameters for Bosonic Excitations**: Fitted parameters of Eq. (5.2) used to model the spectrum in red in Fig. 5.9. n is the number of the bosonic excitation, ω<sup>n</sup> its energy divided by ∆ = (277 ± 4) µeV and a + <sup>n</sup> and a − n are the intensities for positive and negative bias range.


<sup>∗</sup>This mode was manually set to this energy from comparison to the experimental curve. It is found at a too high energy (∼ 1.8∆) by the fit because the BCS fit exceeds the experimental curve.

In the simplest form, the resultant dI/dU spectrum of this excitation ladder can be modeled by fitting a BCS DOS to the gap for E ≲ 2∆, which should constitute the elastic tunnelling signal σel, and then adding this functional form shifted by the respective energy of the inelastic excitation with a variable intensity an. For N inelastic conduction channels involving bosons with energy ωn, the fit function shown as a red line in Fig. 5.9 reads

$$\sigma\_{\text{tot}}(U) = \sum\_{n=0}^{N} a\_n \sigma\_{\text{el}}(U - \omega\_n/e) \tag{5.2}$$

with a<sup>0</sup> = 1 and ω<sup>0</sup> = 0 for the purely elastic part. This approach corresponds to a boson spectrum consisting of infinitely sharp modes, so single δ-distributions at ωn. The experimental curve (black line) was fitted to this model separately for the positive and negative bias range. For U > 0, the bosonic energies ω<sup>n</sup> and intensities a<sup>n</sup> were left as adjustable fit parameters. For the negative bias range, the bosonic energies, which were found for the positive bias range, were enforced, but the intensities were still left adjustable. This choice was totally ambiguous and has also been performed in reverse yielding very similar bosonic energies. As obvious from the measured spectrum, the peaks appear at voltages symmetric with respect to zero bias. Their intensities, however, are different for positive and negative bias voltage. The shape of the spectrum is nicely followed by the fitted function (red line) taking into account N = 5 boson excitations without additional lifetime broadening. The curves in blue, green, violet, yellow and light blue show the fit up to the n'th excitation. The parameters of the fit function are summarized in Table 5.2.

Apart from the second excitation, shown as a violet line, all excitation energies are roughly equidistantly spaced in the spectrum. This is indicated by the dark grey dashed lines in Fig. 5.9. These energies are (within the margin of error) multiples of ωexc ≈ 1.4∆. A dissipative mode of the superconductor at this energy was also found in THz absorption experiments on films of similar sheet resistance, as shown in the inset of Fig. 5.9. The above-mentioned mode is labeled ω ′′. While the modes at ω ′ < ∆ were assigned to two-dimensional plasma phase modes [153, 184], the origin of the resonance at ∆ < ω′′ < 2∆ remains unclear. A phononic origin of the mode at energy ωexc ∼ ω ′′ can be ruled out as van Hove singularities in the phonon DOS of Al or AlO<sup>x</sup> all lie above 15 meV [185]. The same holds for the vibrational modes of chemisorbed oxygen on the aluminium surface [186]. In 3D bulk superconductors the NG mode is usually gapped up to the plasma frequency ω<sup>p</sup> (so to much higher energies than ∆). The drastic reduction to sub-gap energies in grAl samples with high sheet resistance, is ascribed to the granularity [184, 187] and the quantisation of two-dimensional plasmon modes [153] that have maximum DOS at the effective plasma frequency<sup>3</sup> . While bulk plasmonic excitations are expected above ωp, the sharp peaks in energy seen in this work speak more for a non-dispersive mode with a long lifetime.

Recent theoretical studies predict, that the mass of the Higgs mode can be decreased below 2∆ in disordered superconductors [189],making it another plausible candidate. In some of these disordered superconductors near an SIT, Higgs modes with long lifetimes and a renormalized mass just below 2∆ were reported [64, 190]. The mechanism, that was predicted in numerical calculations is, however, still missing confirmation by an analytic theory [191]. It has been shown that in the strong coupling limit, which we expect to be in beyond the Mott resistivity, amplitude fluctuations become equally important as phase fluctuations [98]. The energy scale at which this happens is T ∼ ωexc, which is for the excitation energy determined here, roughly twice the critical temperature. This might be an indicator for the existence of an additional temperature scale Tpair > Tc, in which pairing is already present locally but not in a coherent manner due to large fluctuations. Similar to the case of the cuprates, this would induce a pseudogap above Tc. Such a pseudogap was already proposed as an explanation for the anomalous behaviour of ∆ and the real part of the dynamical conductivity in THz spectroscopy of oxygen-rich grAl samples [149]. Further discussion on the origin of this long-lived mode is postponed to Sec. 5.5 and 5.6, in which our tunnelling data is compared to the predictions within a concrete model for plasmon modes.

In a more general approach, the dI/dU spectrum from Fig. 5.9 was also separated into an elastic and inelastic part and the bosonic spectral density function α <sup>2</sup>F(ω) was determined analagously to how the spin spectrum was obtained for the cuprate superconductors in Chap. 4. This approach has the advantage that the functional

<sup>3</sup> In 2D, plasmons are not gapped [188].

**Figure 5.10: Boson Spectral Density**: Top panel: A separation of the experimentally determined dI/dU spectrum (black line) into an elastic part (red) and an inelastic part (grey) allows for a deconvolution by means of FFT and Gold algorithm. The forward calculated conductance spectrum using the bosonic spectral density function obtained from FFT/Gold algorithm is shown in orange/green. Bottom panel: Bosonic spectral density function obtained via FFT/Gold algorithm (orange/green). Black lines indicate ∆ and 2∆. Dashed lines indicate the energies n · ωexc. Dashed-dotted lines indicate additional resonances at ω <sup>∗</sup> and 2∆ + ω ∗.

form of the bosonic spectrum is not limited to a series of δ-distributions and the number of modes is a result of the deconvolution, not postulated a priori. First, the spectrum from Fig. 5.9 was symmetrized around zero bias in order to enhance the sensitivity for even components of the spectrum, i.e. the particle-hole symmetric parts including the peaks of inelastic excitations. For detailed information on how the data is prepared and the deconvolutions are performed, the reader is referred to Chap. 4 and Appendix A. The symmetrized spectrum is shown in Fig. 5.10 in the top panel (black line) and a BCS fit to the gap (downscaled to an intensity which is always lower than the measured curve) represents the elastic tunnelling conductance (red line). Subtracting this elastic conductance from the experimental curve yields the inelastic tunnelling conductance (grey area) which is a convolution of the bosonic spectral density and the elastic conductance. The determination of the bosonic spectral density α <sup>2</sup>F(ω) was performed via direct fast Fourier transform and the Gold algorithm and is displayed in the bottom panel in orange and green respectively. The respective results of the forward convolution using this bosonic function yields the total conductance in orange/green in the top panel.

Like in Chap. 4, the bosonic spectrum acquired through the FFT deconvolution (orange) contains non-zero values for E < 0, that are cut off in the linear forward convolution. This leads to a calculated conductance below the experimental curve. Also, once again, unphysical negative values are obtained in several energy regions including a region around ∆. Despite these shortcomings, the obtained bosonic spectrum in orange already reveals three broad but prominent modes centred at 0.63 meV, 1.12 meV and 1.51 meV and the experimental dI/dU spectrum is qualitatively well reproduced with it. The bosonic spectrum obtained via Gold algorithm (green) is again limited to strictly positive values in agreement with the physical constraint. The three above-mentioned modes are also found at the same energies, yet with smaller width. Around the peak at 0.63 meV, additional side peaks at 0.455 meV, 0.525 meV and 0.77 meV are found. One low energy mode at 0.07 meV is found at a position where the FFT curve also exhibits a small peak. Although slightly above the experimental conductance curve, the forward convolved dI/dU curve nicely reproduces the observed peaks in position and shape up to 2 meV. A closer look at the mode energies in relation to the gap size ∆ reveals that the peak energies ω<sup>2</sup> = 0.77 meV, ω<sup>3</sup> = 1.12 meV and ω<sup>4</sup> = 1.51 meV correspond to approximately n = 2 , 3 , 4 times ωexc = ω<sup>1</sup> = 0.37 meV = 1.35 ∆, which marks the onset of spectral density between ∆ and 2∆. These energies are therefore connected to n boson excitations with energy ωexc and rising in amplitude with the number n. They are marked by grey dashed lines in the bottom panel of Fig. 5.10. High spectral density at 2∆ is expected to result from efficient pair breaking by inelastic tunnelling processes. The peak amplitude is found at 2∆ + ω <sup>∗</sup> = 0.525 meV, with ω <sup>∗</sup> = 0.07 meV being the sub-gap mode. Since the experiment only yields boson mode energies with no ties to their origin, it remains unclear whether these two modes are connected with each other or not. They are marked as grey dot-dashed lines in the bottom panel of Fig. 5.10. They might also be falsely appointed to the boson spectral density and are really part of the elastic conductivity, i.e. the sample LDOS.

In order to gain more confidence in the presented conclusions on the boson energies and prevent an over-analysis of the data, an independent measurement at a different tip position with different stabilization parameters and bias step size was repeated. The results are shown in Appendix B. The analogous analysis, especially the deconvolution, independently confirms the bosonic mode at ωexc ∼ 1.35 ∆ and even the first excitation is clearly visible there. In total, five prominent excitations are seen, from which the first three fall on multiple energies of ωexc and the last two are slightly shifted to higher energies. Like in the previous spectrum, the weight of the multi-boson excitation seems to grow up to n ∼ 3 − 4 before decreasing again for higher n.

**Figure 5.11: Temperature Dependence of Low Energy Spectrum**: dI/dU spectra at rising temperatures from well below Tc (blue) to above Tc (purple). Individual spectra are offset with respect to each other by 15 nS. Dashed lines indicate the corresponding base line of zero conductance. Data previously published in [146].

At last, the larger pseudogap, which was found with a size of ∼ 8∆ (without fine structure) in a single more insulating grain of sample G1 and with a size of ∼ 14∆ (with a series of equidistant peaks) in the high resistivity sample G2 seems to be shaped by bosonic excitations of the condensate. Whether they are phase or amplitude fluctuations of the superconducting order parameter cannot be concluded from the available data. Such a depletion of states near the Fermi level would happen on a larger energy scale if initiated by the Coulomb blockade in the system. In THz spectroscopy measurements of high-resisitivity grAl samples, the anomalous behaviour of gap size and the real part of the dynamical conductivity has been attributed to a pseudogap above T<sup>c</sup> as a consequence of a phase-fluctuation driven SIT [149]. This is at least in agreement with the results shown here in the sense that this pseudogap only arises in the high-resistivity regime and seems to grow with increasing resistivity. As order parameter fluctuations and long-range order are competing, a coexistance of such a pseudogap and the superconducting gap below T<sup>c</sup> is not too surprising and can be understood as a foretaste of what is to happen at larger grain decoupling at the SIT. How amplitude fluctuations in the small grain regime can destroy superconductivity has been outlined by Guy Deutscher in the 70's [161] and more rigorously explained by Chubukov [98]. Here, measurements above T<sup>c</sup> could be helpful to understand the connection between this pseudogap and superconductivity.

In order to see if inelastic excitations still play a large role above Tc, tunnelling spectra were measured as a function of temperature between 230 mK and 2.19 K. The results are shown in Fig. 5.11 and are offset to each other by 15 nS. The respective base line of zero conductance is indicated by a dashed line for each temperature. Up to 590 mK, the measured grain is superconducting, the gap is fully developed. Above 1.2 K, the measured in-gap conductance does not reach zero anymore due to temperature broadening of the reduced gap. This, however, does not mean that bulk T<sup>c</sup> has already been reached. Between 1.37 K and 2.19 K the spectra exhibit a larger bias range of decreased conductivity that gradually broadens to a V-shaped spectrum at 2.19 K. This is reminiscent of the temperature dependent gap measurements on the cuprates from chapter 4, where the normal state exhibits a similar V-shaped spectrum due to inelastic excitations. Above 2.19 K the spectra do not change significantly anymore which is why the superconducting transition is believed to have happened between 1.37 K and 2.19 K. This is generally in agreement with the transition temperature expected for this sample from the phase diagram in Fig. 5.1. One should keep in mind that this is a local measurement of the gap and variations from heat capacity or bulk transport measurements are expected.

### **5.5** LC **Resonator Model for Plasmon Modes**

That the weight for inelastic tunnelling events including n boson excitations is higher than for the emission of a single boson is surprising but the behaviour of α <sup>2</sup>F(ω)reminds of the multiphoton-assisted tunnelling in superconducting junctions [192]. There, the ratio between microwave amplitude and photon energy yields a cut-off energy that suppresses multiphoton processes with an energy higher than this cut-off and at the same time it irregularly weighs the n multiphoton processes below the cut-off energy in a series of Bessel functions. The emergence of additional side-peaks in dI/dU due to inelastic processes including the absorption and stimulated emission of photons in the microwave field is well described in the Tien-Gordon model [192]. Here, no microwave signal was coupled into the junction. Instead, the logic is reversed: In the absence of a microwave field, only spontaneous emission is possible, which is in agreement with the experimental observation that only high-energy side-peaks are present. The repeated emission of photons in the inelastic tunnelling process leads to energy level fluctuations on the time scale 1/ωexc and with an amplitude Uhf that plays the role of the microwave intensity.

**Figure 5.12:** LC **Resonator Model**: (a) Schematic current flow through the grAl film in the STM. (b) Self-capacitance of each grain is approximated by a metallic sphere with radius R<sup>1</sup> in a dielectric shell of permittivity ε<sup>r</sup> and thickness (R<sup>2</sup> − R1). (c) Electrodynamic model of the network of grains. The tunnelling current It flows through the central line of Josephson junctions in which each grain is capacitively coupled in parallel to other grains that lie in the same plane. (d) 2D plasmon modes for m = 1, 2.., M and different in-plane capacitive couplings C0. Right: Zoom-in on the mode m = 1.

Since the emission of the boson most likely happens inside the sample, one should not think of bare photons like in the vaccum. The sample is a dielectric, thus, the photons strongly couple to the charge degrees of freedom to form plasmons. Before we continue down this road and model the footprint of spontaneous multiplasmon emission processes in the tunnelling spectrum (Sec. 5.6), it is first time to consider a concrete model that legitimates the interpretation as plasmon modes. The main argument, here, will revolve around the energy of the mode.

Fig. 5.12(a) schematically shows how the grAl film is electrically wired in the STM geometry and in which direction the tunnelling current I<sup>t</sup> flows. In our case, the measures of the film are l = 2b = 10 mm and d = 50 nm. We adopt the electrodynamic model from Ref. [187] and describe the grAl film as a network of effective Josephson junctions (JJ). In the RCSJ (resistively and capacitively shunted junction) model [46, pp. 202], each JJ owns a self-capacitance C<sup>J</sup> and a Josephson inductance L<sup>J</sup> . In the Coulomb blockade regime, the aluminium grains act as individual capacitors. Let us assume that each grain is completely embedded in the dielectric with permittivity ε<sup>r</sup> and closely surrounded by other grains which constitute the opposite electrode (Fig. 5.12(b)). We thus model the grain network as N metal spheres of radius R1, each surrounded by a dielectric shell of width (R<sup>2</sup> − R1). The self-capacitance of such a spherical capacitor is then given by

$$C\_J = 4\pi\varepsilon\_0\varepsilon\_r \frac{R\_1R\_2}{R\_1 + R\_2}.\tag{5.3}$$

Evaluating this forR<sup>1</sup> = 1.25 nm,R<sup>2</sup> = 1.75 nmandε<sup>r</sup> = 9.0 [193] yieldsC<sup>J</sup> = 0.73 aF. This value is very close to the capacitance determined through the charging peaks. The Josephson inductance can be determined through the critical current by the formula [30]

$$L\_J = \frac{\Phi\_0}{2\pi I\_c}.\tag{5.4}$$

We use a critical current density of j<sup>c</sup> = 0.88 mA/(µm)<sup>2</sup> for the film with DC resistance ρdc ≈ 2000 µΩ cm [187] and an effective current-carrying area of A = a 2 , where a = R<sup>1</sup> + R<sup>2</sup> is the mean intergrain distance. With these values, we obtain L<sup>J</sup> = 42 nH. The effective plasmon frequency of one JJ,

$$f\_p = \frac{\omega\_p}{2\pi} = \frac{1}{2\pi} \sqrt{\frac{1}{L\_J C\_J}}\,,\tag{5.5}$$

is then f<sup>p</sup> = 914 GHz.

In Fig. 5.12(c), the wiring from Fig. 5.12(a) is translated into an effective circuit model. Each shaded blob depicts a single grain. Since the tunnelling current is injected in z-direction, the current only flows through a handful of grains, i.e. approximately M = d/a effective JJs are put in series. But in the plane of the film, each grain is still capacitively coupled to its nearest neighbours. There are approximately K = b/a grains in a volume V = a 2 b and L = l/a grains in a volume V = a 2 l. We again make use of the electrodynamic model from Ref. [187] and treat the circuit as quasi-1D. This approximation is valid due to the non-uniform electric field induced by the tip electrode. The electric field quickly falls off in the lateral direction of the film, i.e. in first order, the current only flows through the central line in Fig. 5.12(c). The capacitive coupling to neighbouring grains in each plane of the film (vertical connections in (c)), is captured in a mean self-capactiance of each plane C<sup>0</sup> and mainly determined by the average number of nearest neighbours. As demonstrated in Ref. [187], one obtains the energy levels of the quantized 2D plasmon modes by the equation

$$
\omega\_m = \frac{m\pi}{M} \left( L\_J \left( C\_0 + \frac{m^2 \pi^2}{M^2} C\_J \right) \right)^{-1/2} \tag{5.6}
$$

with m = 1, 2, .., M and ω<sup>∞</sup> = ω<sup>p</sup> = 1/ √ LJC<sup>J</sup> being the plasma frequency of the JJ. In Fig. 5.12(d), the energies of the discrete 2D plasmon modes due to confinement in z-direction are shown for C<sup>0</sup> = 2C<sup>J</sup> , 3C<sup>J</sup> , 4C<sup>J</sup> , i.e. two, three or four nearestneighbour grains in each plane. With f<sup>p</sup> = 914 GHz, the bulk plasmons are roughly ten times higher in energy than the bosonic mode we see in our tunnelling spectrum. The first standing wave mode (m = 1) due to confinement by the film thickness d, however, lies between f<sup>1</sup> = 86–121 GHz for C<sup>0</sup> = 2–4C<sup>J</sup> . This corresponds to ω<sup>1</sup> = 355–499 µeV ≈ 1.28–1.80 ∆, which is in the range of our observed ωexc. It should be noted that this is only a rough estimation. The large uncertainty of the critical current density [187] is the major factor in error estimation of ωp, as it only enters L<sup>J</sup> . The grain size plays a secondary role, as L<sup>J</sup> scales inversely but C<sup>J</sup> linearly with the cross-section, such that the dependence nearly cancels out. Additionally, the good agreement between the C<sup>J</sup> estimated from the microscopic model and the value deduced from the charging peaks give reason to believe that the estimate for C<sup>J</sup> is more accurate than for L<sup>J</sup> .

In the following section, the spontaneous multiplasmon emission process is simulated in the Tien-Gordon model. Circumstances under which the excitation of plasmons are easy would be consistent with a low intergrain coupling and low superfluid stiffness, i.e. the system is highly susceptible to phase fluctuations of the order parameter near the SIT [149, 184].

# **5.6 Tien-Gordon Model for Spontaneous Plasmon Emission**

The Tien-Gordon model [192] traditionally describes photon-assisted tunnelling in superconducting junctions that are exposed to a radiation field Uhf cos(Ωt). The electrons can tunnel inelastically via the spontaneous absorption and stimulated emission of photons with the energy Ω. The resultant tunnelling spectra exhibit n sidebands corresponding to the absorption/emission of n photons of energy Ω [194]. How many sidebands occur, depends on the amplitude of the applied field Uhf. Usually, spontaneous emission is neglected because the rate for stimulated emission is much higher in the radiation field.

In this work, no high frequency signal was coupled into the tunnel junction. Instead, the system is considered to be inherently susceptible to inelastic tunnelling events that include the spontaneous emission of n photons in the superconductor. Inside the superconductor, the coupling between charge carriers and photons leads to collective excitations, the plasmons. Via the Anderson-Higgs mechanism, these

**Figure 5.13: Spontaneous Plasmon Emission**: The effective tunnel Eliashberg Function α <sup>2</sup>F(ω) (a) and the differential conductance dI/dU (b) are shown in the analytic description of the Tien-Gordon model considering only plasmon emission. For 5∆ ≲ e|Uhf ≲ 6∆, inelastic processes with up to n = 5 plasmons contribute significantly to the tunnel current and n = 3 is the process with the largest weight (cf. peak height in (a)).

plasmonic excitations also couple to the superconducting phase. It is assumed that the spontaneous emission can, in this case, be described analogously to the spontaneous absorption by simply adjusting the energy conserving term in the δ-distribution. A more rigorous treatment of this problem by the authors of Ref. [195] shows that the total probability for a spontaneous n plasmon emission is described by the convolution of two distribution functions: P0(E) and Pocc(E). P<sup>0</sup> describes the probability that the tunnel junction emits a photon into the vacuum and Pocc describes the absorption capability of the "cavity", which in this case is the grAl film. If the cavity is driven into a coherent state |α⟩ with amplitude α0, then Pocc yields the Tien-Gordon distribution. Additionally, one realizes that the classical AC voltage amplitude Uhf is proportional to the amplitude of the plasmon excitation |α0| [195]. We thus consider our inelastic tunnelling current to behave like a classical drive to the resonating aluminium grains.

Neglecting absorption, i.e. considering only the high-energy sidebands, and restricting ourselves to a single plasmon mode with energy ωexc, the total conductance for single quasiparticle tunnelling in the Tien-Gordon model is written as

$$\sigma\_s^{\rm tot}(eU) = \int\_0^{eU} \mathrm{d}\omega \sigma\_s^{\rm el}(eU - \omega) \sum\_n J\_n^2 \left(\frac{eU\_{\rm hf}}{\omega\_{\rm exc}}\right) \delta\_a(\omega - n\omega\_{\rm exxc})$$

$$= \sigma\_s^{\rm el}(eU) + \int\_0^{eU} \mathrm{d}\omega \sigma\_s^{\rm el}(eU - \omega) \underbrace{\sum\_{n=1} J\_n^2 \left(\frac{eU\_{\rm hf}}{\omega\_{\rm exc}}\right) \delta\_a(\omega - n\omega\_{\rm exc})}\_{\propto \alpha^2 F(\omega)}\tag{5.7}$$

with

$$\delta\_a = \frac{1}{|a|\sqrt{\pi}} \mathbf{e}^{-(x/a)^2}. \tag{5.8}$$

J<sup>n</sup> are Bessel functions of order n and Uhf is proportional to the amplitude of the plasmon state. The Bessel function for n = 0 describes elastic tunnelling. Hence, the terms with n > 0 can be collected to an effective "tunnel Eliashberg function" α <sup>2</sup>F(ω).

In Fig. 5.13, the function α <sup>2</sup>F(ω) and the resulting differential conductance in the Tien-Gordon model is shown for ωexc = 1.35∆, σ el s (eU) = σ0|eU|/ p (eU) <sup>2</sup> − ∆<sup>2</sup>, a = 0.2∆ and the three indicated values for Uhf. All α <sup>2</sup>F(ω) were normalized by the condition of the elastic conductance

$$
\sigma\_{\rm s}^{\rm el}(eU) = J\_0^2 \left(\frac{eU\_{\rm hf}}{\omega\_{\rm exc}}\right) \int\_0^{U} \mathrm{d}\omega \sigma\_{\rm s}^{\rm el}(eU-\omega)\delta\_a(\omega) \stackrel{eU \to \infty}{\rightarrow} \sigma\_0. \tag{5.9}
$$

The ratio eUhf/ωexc yields a cut-off for the number of plasmons that are excited in one inelastic tunnelling event. For eUhf ≲ 6∆ the maximum number is n = 5. Due to the functional form of the Bessel functions, the ratio eUhf/ωexc also adjusts the relative amplitudes for the multi-plasmon excitations.

As can be seen from Fig. 5.13, good agreement with the experimental data, i.e. the process n = 3 has the largest weight, is found for 5∆ ≲ eUhf ≲ 6∆. The resulting dI/dU spectrum in the Tien-Gordon model then closely resembles the experimental spectrum with the characteristic excitation ladder up to a cut-off at e|U| ≈ ∆+5·ωexc.

## **5.7 First Results on Oxidized Polycrystalline Al**

The STS measurements of the YSR states in oxygen-rich grAl (sample G1) from section 5.3 suggested that it is likely that the unpaired spins reside in the nonstoichiometric oxide of aluminium. In order to test this hypothesis, the pure and surface clean aluminium film (sample A) from section 5.2.1 was oxidized by exposure to a pure O<sup>2</sup> atmosphere. Oxygen gas (purity= 99.9999 %) was introduced to the UHV chamber through a variable leak valve up to a partial pressure of p<sup>O</sup><sup>2</sup> = 1 × 10<sup>−</sup><sup>5</sup> mbar and th sample was oxidized for 1 h at room temperature. This corresponds to a dose of roughly 2.7 × 10<sup>4</sup> L (Langmuir). Due to the quick passivation of the aluminium surface after the initial incorporation mechanism (the first chemisorbed monolayer of O<sup>2</sup> is buried in the second surface layer by updiffusion of aluminium), the oxygen uptake rises only slowly for higher doses [196] and saturates after roughly 2 hours [197]. The exposure dose used here should lead to no more than 2 ML of AlO [197–199]. If dangling bonds in grain boundaries and interfaces of the oxide are the host for unpaired electrons, then YSR states should also be observed at the surface of the oxidized aluminium film. Conveniently, complicated factors like size effects of the Al grains, disorder of the superconductor and percolation effects are taken out of the equation in this control experiment.

**Figure 5.14: Oxidized Al Film**: dI/dU spectrum with superconducting gap obtained on the strongly oxidized Al film at T = 39 mK. No YSR states are visible. The blue line shows a Dynes DOS fit with ∆ = (128 ± 1) µeV and quasiparticle scattering rate Γ = (65 ± 1) µeV. Feedback condition: U = 0.5 mV, I = 100 pA. Inset: Large bias range I(U) curve displaying the low conductance of the junction. Due to the large insulating gap the STM tip is in contact with the surface for low bias voltages. The tunnel barrier is the aluminium oxide not the vacuum.

Fig. 5.14 shows the dI/dU spectrum of the oxidized aluminium film at T = 39 mK. The superconducting gap is reduced and not fully developed. The coherence peaks are also broadened significantly. In-gap states are not visible. The modified appearance of the superconducting gap is caused by the fact that tunnelling into the superconducting aluminium is only made possible by bringing the STM tip in contact with the surface oxide. In contrast to the vacuum, the oxide layer allows not only ballistic but also diffusive transport. That means that the tunnelling electron can in principle inelastically scatter on defects inside the oxide. This form of dissipation is seen as excess sub-gap differential conductance and absent coherence peaks. The exponential decay of the tunnelling current makes it virtually impossible to record a tunnelling current if the tip is free standing above the surface because an additional insulating layer of ≈ 0.8 nm is separating it from the superconducting aluminium. The good insulation of the thin surface oxide is well captured in the large range I(U) curve in the inset of Fig. 5.14. Only at a stabilization voltage of U > 2 V a tunnelling current of few pA starts to flow. To account for scattering processes of electrons with sub-gap energies in the presence of the oxide, the superconducting gap was fitted to a Dynes DOS function (see Chap. 4 for a functional expression). The fit yields ∆ = (128 ± 1) µeV and a scattering rate of Γ = (65 ± 1) µeV. If YSR states are not present or just overshadowed by the poisoning scattering channels for sub-gap energies cannot be discerned.

This seemingly simple control experiment beared unexpected complications due to the quick evolution of a large insulating gap of roughly 4 eV, already for a supposedly thin oxide layer. A great amount of preparation of a tip with good energy resolution (in order to resolve the YSR states) is opposed to its one time use. As soon as the tip is dipped into the oxide, the measurement has to be performed at once. At each lift off the surface the tip picks up a cluster of the insulating oxide, which successively increases the effective tunnelling barrier for the next spectroscopy until no tunnelling current can be detected anymore. It is expected that the measurement of an oxidized Al film with a thickness of 0.5 ML < Θ < 1 ML yields more conclusive results concerning the existence of YSR states at the AlO surface. To study the depencence between oxide thickness and abundance of these states, however, the experimental procedure would have to be adjusted to eliminate the problems encountered in the tunnelling geometry chosen here.

Density functional theory (DFT) calculations of chemisorbed O<sup>2</sup> on a few layers of α−Al2O<sup>3</sup> on Al(111) predict that the charge transfer of the oxygen molecule heavily depends on the oxide thickness and an effective spin of 1 µ<sup>B</sup> on the adsorbant is expected only at 4−9 trilayers of Al2O<sup>3</sup> [200]. In order to realize such thick aluminium oxide layers either much larger oxygen pressures [199] or elevated temperatures that facilitate the surface reactivity (dissociative chemisorption of O<sup>2</sup> is thermally activated) [201] are needed. The latter method would be the preferred in a UHV environment like the setup used in this work.

## **5.8 Summary**

The combined STM/STS measurements on grAl in the low resistivity and high resistivity regime confirmed the existence of nano-sized grains and yielded enhanced, uniform superconducting gaps of ∆ = (298±1) µeV and ∆ = (312±1) µeV compared to ∆ = (191 ± 1) µeV for a polycrystalline aluminium sample. The oxygen-poor grAl sample displayed metallic behaviour for larger bias voltages, which speaks for well coupled grains, and the superconducting gap was free of sub-gap states. In the oxygen-rich grAl sample, charging effects, indicative of the Coulomb blockade regime, were observed. Apparently, the grains decouple as a function of oxide thickness and can be collectively charged in the form of clusters containing several electrically well connected grains. In some grains, YSR states were found as a fingerprint of localized spins. The energetic shift of the YSR states in dependence of an applied electric field suggests that the unpaired spins reside in the oxide barriers. An independent STS experiment on a polycrystalline aluminium film with an estimated surface oxide thickness of ≈ 2 ML revealed no YSR states. It is likely that the localization of unpaired spins requires a thicker oxide. On individual grains in the oxygen-rich grAl sample, the out-gap structure of scanning tunnelling spectra is shaped by inelastic excitations with long lifetimes. The extracted effective bosonic spectral density of these excitations can be qualitatively well explained by the spontaneous emission of multiple plasmons. This would very much support an SIT that is driven by phase fluctuations of the superconducting gap. From an energy standpoint, we argued that the excited plasmon mode is most likely part of the lowest branch of discrete modes that are quantized through the thickness of the grAl film. An alternative, yet very speculative, candidate for these inelastic excitations are Higgs modes.

# **6 Multiflux Vortex States inBulk Pb**

The recent evidence for two-band superconductivity [39] and its naturally close proximity to the Bogomolnyi point [41, p. 436] make lead a prime candidate to study the interplay of multiband superconductivity and normal conducting domains in the intermediate state. Bulk multiband superconductivity in conjunction with conventional s-wave pairing has been proposed to unveil a new route to topological superconductors that are robust to magnetic fields and host Majorana modes [202, 203]. At medium magnetic fields, the low-temperature transitional region of a multiband superconductor with close to critical κ, is also predicted to host interesting structures like multiflux vortices [204]. These vortices represent topological objects themselves, that may also host Majorana modes under the right conditions [205]. In this chapter, new findings on the appearance of the low-temperature intermediate state of Pb are presented. These include the emergence of single-flux vortices and multiflux vortices. In the process, the influence of interband coupling, tunnelling into Caroli-de Gennes-Matricon states and a robust method for winding number determination are detailed. The numerical simulation of the vortex patterns via solution of 3D Eilenberger equations confirms the interpretation of presented scanning tunnelling data. The Fermi velocities that enter the numerical simulations were provided by Rolf Heid from IQMT, KIT, and are obtained within a fully relativistic DFT approach [206]. This chapter is based on the following preprint: [207].

### **6.1 Meißner and Intermediate State**

After several cycles of hot sputtering and annealing, as described in Sec. 3.3.2, only very few adsorbants remained at some kink sites of step edges on the Pb(111) surface and terraces formed with widths of up to few hundred nanometres, as can be seen in the topographic scan shown in Fig. 6.1(a). Upon zero-field cooling the Pb(111) sample to 39 mK, it enters its superconducting state below T<sup>c</sup> ∼ 7.2 K [41, p. 436]. The dI/dU spectrum in Fig. 6.1(b) demonstrates that the two superconducting gaps of Pb can be resolved in our setup, even with a normal conducting tip. From the fit of a temperature broadened two-gap BCS function (red curve), the gaps are

**Figure 6.1: Superconducting State**: (a) Topographic scan of the Pb(111) surface showing almost no contamination and wide terraces. (b) Differential conductance spectrum at T = 39 mK centred around the Fermi energy. The superconducting gaps ∆1,<sup>2</sup> are determined from a broadened two-gap BCS function fit (orange curve). The insets show 3D models of the FSs of the two superconducting bands (from Ref. [209]) inside the first BZ (wireframe) viewed from the (111) direction. ∆<sup>1</sup> is assigned to the band forming a tubular FS (red). The larger gap ∆<sup>2</sup> is assigned to the band with a compact FS (cyan).

determined as ∆<sup>1</sup> = (1.26 ± 0.02) meV and ∆<sup>2</sup> = (1.40 ± 0.02) meV. This is in good agreement with previous measurements of the difference of the two gaps [39], where the authors made use of the energy resolution amplification of superconducting tips.

It should be noted that the intensity difference of the two coherence peaks has previously been attributed to the k-dependent tunnelling matrix element and the smaller (larger) gap has been assigned to the compact (tubular) Fermi surface (FS) sheet [39], shown in the cyan (red) 3D model in Fig. 6.1(b). This is in contrast to BdG calculations [208] and the findings in this work. A conclusion from the analysis in Sec. 6.2.2 is that the band-to-gap assignment should be reverse: The band with the tubular/compact FS (red/cyan) is from here on denoted as band I/II and is responsible for the superconducting gap ∆1/∆2.

After applying a perpendicular magnetic field of B = 85 mT, which is above the critical field of µ0H<sup>c</sup> ∼ 80 mT [210], magnetic flux enters the sample from the sides and destroys superconductivity completely. Upon subsequent decrease of the field below Hc, the Landau intermediate state is reached. This state is characterized by large normal and superconducting domains and by long domain walls separating the two. The shapes and sizes of domains in the intermediate state of lead have been extensively studied by magneto-optical methods revealing the strong dependencies on temperature, sample shape and magnetic protocol [211–214]. Since the domain sizes are in the range of several micro- or even millimetres, STM is not exactly suited to study these domains. Still, it is possible to frequently freeze domain walls in the scan frame by performing a slow magnetic field ramp (5 mT/min) and simultaneously measuring the differential conductance at eU<sup>t</sup> ∼ ∆<sup>2</sup> at the tip position. As soon as a

**Figure 6.2: Intermediate State**: (a) dI/dU map at U<sup>t</sup> = 1.3 meV showing a typical domain wall in the intermediate state at B = 23 mT. The normal (superconducting) domain exhibits a low (high) conductance at this bias voltage. Right overlay: Corresponding topographic scan image of the area. (b) Schematic cross-sectional view through the bulk material in the intermediate Landau state. The magnetic flux density (blue) is favouring a more homogeneous distribution at the surface than in the bulk leading to branching. (c) dI/dU spectra along a cross-section from normal conducting to superconducting domain indicated by the white frame in (a). The direction of the profile is indicated by the red arrow. The spectra are locally averaged over a straight part of the domain boundary (perpendicular to the profile direction) and recorded in distance increments of ∆d = 9.19 nm.

domain wall passes the tip position, the measured conductance jumps and the ramp is stopped. After a few seconds of waiting until the dynamics have stopped, finding a domain wall within a scan frame of 1.4 µm × 1.4 µm is likely.

A typical domain wall in the intermediate state is shown in the dI/dU map in Fig. 6.2(a). At a tunnelling bias of U<sup>t</sup> = 1.3 mV, the normal conducting domains show up as areas of low conductance (purple) and the superconducting domains as areas of high conductance (green/yellow). Thin lines of high conductance are an artefact of the measurement mode and appear at the position of step edges on the surface as the right overlay shows. They are the consequence of a sudden jump in the tunnelling current when the tip gets too close to an edge in the multipass mode. They are therefore ignored in the discussion of the LDOS. A cross-sectional line scan across the domain wall, as in Fig. 6.2(c), shows how the gap ∆<sup>1</sup> changes from zero to its maximum on the length scale of the coherence length. The local recovery of superconductivity agrees well with reported coherence lengths of ξ ∼ 87 nm [41, p. 436]. A detailed analysis of the coherence length ξ1,<sup>2</sup> of the two bands, measured inside vortices, will follow in Sec. 6.2.1. Direct information about the magnetic field lines and therefore also the local supercurrents inside the superconductor bulk cannot be gained from neither STM nor the magneto-optical (MO) method mentioned earlier. From energetic considerations, however, the intermediate state is supposed to form according to two competing effects: The interface energy to build domain walls is minimized by a course distribution of magnetic flux in the interior, while the magnetic field outside the sample favours homogeneity. The resulting compromise is a branching flux pattern, as shown in Fig. 6.2(b), which was first predicted by Lev Landau [215, 216]. These thermodynamic considerations neglect the influence of flux pinning centres in real crystals and dynamic effects caused by the magnetic field ramp, Lorentz force acting on the flux tubes and topological hysteresis [212].

### **6.2 Emergence of Isolated Vortices**

For disk-shaped Pb samples, R. Prozorov could show that the equilibrium state after flux penetration of a zero-field cooled (ZFC) sample is tubular: With increasing field, flux enters the sample from the outer edges and is quickly pushed to the centre by the Lorentz force. The tubular structure grows with increasing field in a coarsening froth pattern. At the same field values, the sample exhibits a laminar magnetic flux pattern on the surface (in accordance with Landau's prediction) if this state has been reached from the normal state upon field decrease [212]. While this is an experimental observation on the micrometer scale and at 5 K, it is nevertheless highly unexpected that isolated vortices, so flux tubes carrying a single flux quantum, form in Pb at ∼ 30 mK when the field is decreased following the magnetic protocol described in the previous section. If single-flux vortices are indeed stable in a bulk Pb single crystal, then it raises the question whether strong pinning centres can explain this fact, or the classification of superconductors into type-I and type-II has to be revisited at ultra-low temperatures T ≪ Tc.

#### **6.2.1 Normal Single-Flux Vortices**

Here, compelling evidence for the existence of stable magnetic vortices in intermediate state bulk Pb(111) is shown by display of the local density of sub-gap states in a magnetic vortex, called Caroli-de Gennes-Matricon (CdGM) states. The quantum mechanical nature of theses states is explained in detail in Sec. 6.3. For now, the following quasi-classical picture suffices to discuss most features of the experimentally observable LDOS: In a magnetic vortex core, the superconducting order parameter or gap size vanishes. Normal state electrons and holes with E < ∆ are trapped in this vortex as they are (mostly) specularly reflected by the potential barrier ∆. Thus, a great number of quantum well states (QWS) exist inside the magnetic vortex core. Their abundance at a certain energy is determined by the normal state DOS and their penetration depth into the superconducting surrounding is determined by the Fermi velocity v<sup>F</sup> 1 .

After ramping the magnetic field down from 85 mT, isolated, round normal conducting domains of only ∼ 100 nm in diameter are frequently stabilized. They appear as a conductance depression at eU<sup>t</sup> = ∆2, as displayed in Fig. 6.3(h). The appertaining zero bias conductance map is shown in Fig. 6.3(a). It reveals a threefold symmetric star-shaped pattern with a maximum in the star's centre. The CdGM LDOS stretches over 100 nm in the ⟨2¯1¯1⟩ directions in the form of thin rays (blue/green) that continuously fade at increasing distances from the centre. On top of that, three weak star arms (dark blue) along the ⟨1¯21⟩ directions are visible. The bulk crystal directions have been determined from larger crystal defects that reveal the direction of glide planes. Further information can be found in Chap. 7.

With increasing energy (independent of sign of bias voltage), the star's weak and strong arms both split into two with increasing splitting distance (see Fig. 6.3(b-f)). Meanwhile, the central peak splits isotropically to a ring shape (Fig. 6.3(b-f)). At eV ≲ ∆<sup>1</sup> the ring size is maximal (see Fig. 6.3(g)). This pattern proved to be invariant under reversal of tunnelling current or magnetic field direction, which is demonstrated in Appendix D and showcases the particle-hole symmetry of the CdGM states.

By solving the quasi-classical 3D Eilenberger equations, the experimentally found LDOS patterns of CdGM states could be qualitatively reproduced for bands I and II separately. In Fig. 6.3(a,e), the right bottom (top) panel shows the calculated LDOS signature for a single-flux vortex using the DFT calculated band structure from band I (II) at given energy. Apparently, the long, stretched-out star arms of high LDOS stem from electrons at the Fermi edge of band II, for which the FS is shown in cyan. A view onto the FS from the (111) direction makes it obvious, that these long arms are caused by the large Fermi velocity anisotropy in this band. A large portion of the FS consists of flat parts where the Fermi velocity's in-plane component points in one

<sup>1</sup> For a non-relativistic electron, the classical momentum p = mv can be used in the kinetic energy term of the Hamiltonian Hkin(r) = −iℏv∇.

**Figure 6.3: Single Flux Vortex**: Multi-pass images of the differential conductance at sub-gap energies (a-h), recorded at T = 37 mK and Bext = 0 mT. The inset in (a) shows the atomic lattice and important bulk crystal directions are marked by yellow/red arrows. Right panels in (a,e) complement experimental data with LDOS simulations using a 3D Eilenberger formalism and DFT calculated Fermi velocities. The LDOS pattern was calculated seperately for band I (top) and band II (bottom) with FSs as shown in the small inset.

of the ⟨2¯1¯1⟩ directions. The splitting of the arms is also seen in the semi-classical calculation at non-zero energy. The root of this effect lies in the quantum-mechanical nature of these states and will be discussed in detail in Sec. 6.3. The calculations for band I confirm that the weak arms in the ⟨1¯21⟩ directions as well as the strong central LDOS peak stem from electrons at the red FS. The openness of this FS provides for a larger variety of Fermi velocity direction, which effectively leads to a more isotropic scenario, albeit the threefold symmetry still remains apparent. The stronger confinement of states in the centre is a result of the (on average over the whole FS) larger z-component<sup>2</sup> of the Fermi velocity. Consequently, the respective

<sup>2</sup> Here, the z-direction denotes the direction perpendicular to the crystal surface, so the [111] direction.

in-plane component is smaller, which leads to a smaller penetration depth into the superconducting surrounding. The out-of plane component of the Fermi velocity proved to play a crucial role in the Eilenberger calculations, as solving effective 2D Eilenberger equations in the (111) plane would only yield sixfold symmetric LDOS patterns. The 3D trajectories of the quasiparticles do, however, have substantial z-components. Upon respecting the translation symmetry breaking by the surface, in conjunction with 3D quasiparticle trajectories, the threefold symmetric LDOS patterns could finally be reproduced. A detailed explanation of the simulation methods used, can be found in Appendix C.

It should be noted that this is not the first observation of anisotropic LDOS signatures due to CdGM states in and around a vortex structure: The first to observe such star-shaped structures were Hess et al., who observed them in the vortex lattice of 2H-NbSe<sup>2</sup> [217]. Later, Hayashi et al. demonstrated that the semi-classical theory of Eilenberger can reproduce the experimental results nicely [218, 219]. Recently, they were also observed on the (0001) surface of the elemental superconductor La [220], where the anisotropy was attributed to the Fermi velocity and not to an anisotropy in the gap parameter which was shown to have a large affect the LDOS pattern, too [221, 222]. This work marks the first experimental observation of highly anisotropic CdGM states from two different bands inside solitary vortices. The early stated presumption that the anisotropic pattern might be caused by neighbouring vortices [217, 223] in a flux lattice can therefore be ruled out. In addition, this is the first time, a pattern with odd-integer rotational symmetry is found. It could be shown that this pattern can be reproduced using 3D Eilenberger simulations with accurate Fermi velocites and that the total signature can be understood as a superposition of the CdGM states from both superconducting bands.

After learning which sub-gap states belong to which band, it makes sense to take a deeper look into the single bias spectroscopies inside the vortex in order to determine the length scale on which ∆<sup>1</sup> and ∆<sup>2</sup> change. Previously, in Sec. 6.1, this proved to be difficult for the domain wall due to the ripples inside it<sup>3</sup> . Radial linegrids recording the differential conductance were taken on the vortex from Fig. 6.3. In Fig. 6.4, the second derivative of the differential conductance −∂ <sup>2</sup>σ/∂V <sup>2</sup> is shown for four different directions (indicated in the inset). The grey color code has been adjusted such that white regions mark clear peaks in the differential conductance curves. One can now follow the traces of these peaks as the distance from the centre increases. For α = 0◦ the linegrid is aligned with a strong star arm. We can thus recognise the trace of peaks at zero energy to be the CdGM states of band II (cyan

<sup>3</sup> The necessary averaging over a certain length of the domain wall had a negative impact on the overall energy resolution and made a distinction between ∆<sup>1</sup> and ∆<sup>2</sup> difficult over a large distance.

**Figure 6.4: Coherence Lengths in the Normal Single-flux Vortex**: Angle dependent radial measurement of the vortex core states from the vortex in Fig. 6.3. Displayed is the second derivative of the differential conductance −∂ <sup>2</sup>σ/∂V <sup>2</sup> in order to highlight maxima in the LDOS. The point distance of single spectra is ∆d = 2.57 nm. Marked are the opening of ∆<sup>1</sup> (red line), the opening of ∆<sup>2</sup> (cyan line) and the CdGM states of band II (cyan dashed line). The inset shows the direction of the y-axis.

dashed line). For 0 ◦ < α < 60◦ the states move to higher energies and eventually fall together with the opening of ∆<sup>2</sup> at α = 60◦ (cyan line). A crossing of ∆1(r) and ∆2(r) is identified at d ≈ 50 nm. That means, that the length scales on which the two recover to their maximum value must be different. The smaller gap ∆<sup>1</sup> seems to recover much quicker than ∆<sup>2</sup> when moving away from the vortex centre (red line). While ∆2(r) can be described by ∆2(r) ∼ ∆<sup>2</sup> tanh(r/ξ2) with ξ<sup>2</sup> ∼ 45 nm over the whole distance, ∆<sup>1</sup> cannot. It follows a tanh behaviour with ξ<sup>1</sup> ∼ 45 nm for large distances from the centre but has a much steeper slope yielding a core size of only ξ (c) <sup>1</sup> = ∆1(∞) h limr→<sup>0</sup> d∆(r) dr i<sup>−</sup><sup>1</sup> ≈ 10 nm.

Such a shrinking of the vortex core size is known from theory as the Kramer-Pesch effect [223, 224]. It is so far mainly studied in isotropic vortices [224–226] and the self-consistent solution of the pair potential with isotropic Fermi velocity shown in Fig. D.2 confirms that a Kramer-Pesch effect should be present and lead to a significant deviation from a tanh shape of ∆(r). This observation is in accordance with our band-gap assignment: Band I, which is responsible for ∆1, has a more isotropic Fermi velocity distribution than band II, for which the Kramer-Pesch effect is missing. It makes sense that, in the anisotropic case, the thermal population of the CdGM states at T ≪ T<sup>c</sup> does not lead to a core size reduction because the low energy CdGM states can extend far from the vortex centre in certain directions. For a superconductor in the clean limit, the core size should shrink proportional to T /T<sup>c</sup> and should even be able to become infinitely small at T → 0 K [227]. This raises the question, how the opening of the gap should rightfully be determined in the low temperature experiment. As one can see in Fig. 6.4, the quasiparticle coherence peak at ∆<sup>1</sup> is already present at ≈ 10 nm distance from the centre which could mean that the core size is even smaller and the red line follows the CdGM states of ∆<sup>1</sup> which do not coincide with the "opening" of ∆1. This is, however, almost a philosophical question: Since these CdGM states are always present in a vortex, they always take away from the condensate and only when they reach E = ∆ is the corresponding coherence peak maximally large and the gap fully empty. Apparently, there is no gap size matching (the similar ∆ are rather a coincidence) and no length scale locking present in Pb. This leads to the conclusion that the bands are relatively decoupled from each other [87]. This readily justifies the separate treatment of the two bands in the quasi-classical simulation.

#### **6.2.2 Anomalous Vortices**

Apart from normal single-flux vortices, anomalous vortices like the one in Fig. 6.5 are frequently observed. These vortices also carry one single flux quantum and behave like the normal vortices but the sets of CdGM states from the two bands seem to be spatially displaced with respect to each other. The fact that the strong central zero bias peak in Fig. 6.3(a) belongs to a different band than the strong star arms was previously only inferred from the comparison with the theoretical calculations but can now be confirmed experimentally. In the anomalous vortex, this strong zero bias peak does not lie in the centre of the star pattern, as can be seen in Fig. 6.5(a). In Fig. 6.5(b-g), it is also much clearer that this peak splits almost isotropically to a ring-like pattern for higher energies. Due to the displacement of the CdGM states, the zero energy LDOS in the star centre and the ring centre can be compared. This is shown in Fig. 6.5(i). The amplitude of the tunnelling conductance in the ring centre is roughly four times larger than in the star centre. In fact, the band-gap assignment is unambiguously possible by using the fact that the CdGM states from different bands are spatially displaced: By recording the differential conductance in fine linegrids (Fig. 6.5(j,k)) along the paths k,j, that are marked in Fig. 6.5(i)'s right panel, one sees that ∆<sup>2</sup> (cyan line) closes in the star centre, while ∆<sup>1</sup> (red line) closes in the centre of the ring. Since the CdGM state-band assignment is known from the quasi-classical simulations, one can now also definitively relate band I (tubular FS) to the smaller ∆<sup>1</sup> and band II (compact FS) to the larger ∆2.

The direction and magnitude of spatial displacement between the two state sets could not be tied to any crystal direction and appears to be almost random. By changing the magnetic field, the displacement can be modified, as Fig. 6.6 demonstrates. In fact, both CdGM state sets are mobile and as a result, the anomalous vortex can be transformed back into a normal vortex. As a point of proof, Fig. 6.7(a,c) shows the respective topographic images to the conductance maps in (b,d). Clearly, the scan

**Figure 6.5: Anomalous Single-flux Vortex**: (a-h) dI/dU maps like in Fig. 6.3 for an anomalous vortex. (i) Right panel: Enlarged image of (a) including the position of line grid spectra (white arrows) and single spectra locations (red/black circle). Left Panel: Single dI/dU spectra in the star centre (black) and ring centre (red) revealing zero bias peaks of different amplitude. (j,k) Heat maps of the second derivative of the differential conductance −∂ <sup>2</sup>σ/∂V <sup>2</sup> like in Fig. 6.4 along the cross-sections marked in (i). The red/cyan line follows ∆1/∆2.

area remains the same but the vortex moved and transformed from the anomalous type to the normal type. The chances of this being two entirely different vortices are rather small considering the low probability to find one in the first place. So overall, the CdGM states of individual bands seem to be almost disconnected from each other inside one vortex, which ties into the previous indications of a small interband coupling. Additionally, the vortices themselves are moveable by a change in magnetic field and can be changed from the anomalous to the normal type. One explanation for this could be, that the flux tubes are tilted with respect to the surface and that this tilting causes the displacement of the CdGM state sets by removal of the cylinder symmetry. In fact, the bright stripe (which also splits in energy) that points toward the vortex in Fig. 6.5, 6.6 and 6.7 is only present for anomalous vortices

**Figure 6.6: Manipulation of CdGM State Displacements**: Upon slowly varying the magnetic field, the relative displacement of the star centre and ring centre (red) changes in distance and direction untied to any crystal direction. (a-d) Differential conductance at zero bias voltage around the vortex from Fig. 6.5 at decreasing magnetic fields. The position of the ring centre is marked by a white-dashed circle. The bright stripe mentioned in the main text is marked at its borders with white dashed lines

and might indicate the direction in which the magnetic field lines are tilted. For clarity, it is marked in Fig. 6.6 by white dashed lines, together with the position of the ring centre by a white dashed circle.

**Figure 6.7: Manipulation from Anomalous to Normal Vortex**: (a,c) Topographic images of the same area on the surface before (a) and after (c) the decrease in magnetic field. (b,d) Corresponding dI/dU maps to the topographic images (a,c). After a slow decrease of the magnetic field, the anomalous vortex in (b) moved and transformed into a normal type (d). In the process, both sets of CdGM states moved.

### **6.3 Tunnelling into Caroli-de Gennes-Matricon States**

By solving the BdG equations for a cylinder symmetric vortex line in a type-II superconductor, C. Caroli, P.G. de Gennes and J. Matricon could show that there exist fermionic bound states with energies smaller than the gap size ∆<sup>∞</sup> in the vortex core [229]. For ϵ ≪ ∆∞, these Caroli-de Gennes-Matricon (CdGM) states approximately follow the linear dispersion

$$
\epsilon(\mu, k\_z = 0) \propto \mu \frac{\Delta\_{\infty}}{k\_F \xi} \approx \mu \frac{\Delta\_{\infty}^2}{\epsilon\_F} \tag{6.1}
$$

with µ being an angular momentum quantum number that defines each state |jµ⟩ with total angular momentum j. One could think of these states as quantised Landau levels with the difference being that each of these states has to be built from the superconducting condensate and is therefore a mixture of electron and hole character. The spectrum for a single-flux vortex, calculated by Gygi and Schlüter, is shown in Fig. 6.8(a). The states are in principle discrete in energy, a zero-energy state is missing and the spectrum is inversion symmetric to the origin ϵ = µ = 0, meaning

**Figure 6.8: CdGM states**: Quantum mechanical calculation of ϵ(µ) for vortices of different vorticity m. (a) Energy spectrum for an m = 1 vortex. Reprinted figure with permission from [223]. Copyright (1991) by the American Physical Society. (b,c) Energy spectrum for an m = 4 (b) and m = 5 vortex (c). Reprinted figure with permission from [228]. Copyright (1999) by the American Physical Society.

**Figure 6.9: Multiquantum Vortex Fingerprint**: Calculated radial conductance profile for an m-quanta vortex. Reprinted figure with permission from [228]. Copyright (1999) by the American Physical Society.

that the excitations for µ < 0 have ϵ(µ) < 0. The dispersion is linear for low energy and flattens as the CdGM states approach the gap energy. While individual CdGM states can be probed on a few unconventional superconductors like FeTe0.55Se0.<sup>45</sup> [230], that have a relatively high gap energy and small Fermi energy, for conventional superconductors the ratio ∆2/ϵ<sup>F</sup> is, even at ultra-cold temperatures, below the resolution limit of current techniques like STS. For Pb, this energy spacing is of the order of ∼ 0.1 µeV. As a result, the states form a quasi-continuum, which is called a CdGM branch. As can be seen from comparison with results for an m = 4 and m = 5 vortex (see Fig. 6.8(b,c)) by Virtanen and Salomaa [228], there exist m individual CdGM branches in an m-quanta vortex. While all spectra are inversion symmetric to the origin, only vortices with odd winding number have a branch that crosses E = 0 at µ = 0.

An important feature of these states, that emphasizes their similarity to Landau levels, is that the quantisation translates to some extent into real space: Gygi and Schlüter could show that a CdGM state with angular momentum µ has maximum amplitude at a distance r ∼ |µ|/k<sup>F</sup> from the vortex centre [223]. That is for a cylinder symmetric vortex and an isotropic in-plane Fermi velocity. This makes STM an excellent choice of instrument to study these states, as the whole spectrum of CdGM states can be probed in radial STS profiles of the vortex. These profiles should then exhibit the wedge shape from Fig. 6.9 with multiple dispersing branches for an m-quanta vortex. If all CdGM state branches can be disentangled and identified in these profiles, then an unambiguous determination of the vortex winding number is possible.

#### **6.3.1 Determination of Vortex Winding Numbers**

In the same year, in which Gygi and Schlüter published their results from the microscopic theory, G.E. Volovik made the connection of the low energy fermionic states inside a vortex, that he studied primarily in rotating superfluid <sup>3</sup>He-B, to topology. He found, that the number of zero energy crossings Nzm, i.e. points in k-space where particle and hole spectrum meet, is a topological invariant and therefore very robust to perturbation of the system [227, 231, 232]. This means, that the overall shape of the vortex should not matter and that the number of individual CdGM state branches is simply determined by the vortex winding number m, in particular

$$N\_{\rm zm} = m.\tag{6.2}$$

As discussed beforehand, in the case of the cylindrically symmetric vortex with isotropic Fermi velocity, m can easily be determined in an STM experiment. That is because the number of zero-energy peaks in a radial STS profile going from −ξ to +ξ through the vortex centre is Nzm and thus identical to m. The problem in deviation from the model case of perfectly cylindrical vortices and isotropic bands is that the description of the axial degree of freedom by an angular momentum is not correct, and the solutions of the problem will deviate from simple Bessel functions. Consequently, the quasiparticle spectrum, whether described fully quantum mechanically or semi-classically, has to be described by the quantum numbers (kx, ky, kz). The application of the topological index theorem then proves to be more difficult from an experimental point of view, as it is unknown which line in real space crosses all diabolic points (the zero energy points in the spectrum, in this case) in k-space exactly once. Even though the states do not possess a continuous rotation symmetry anymore, the next section will show that the shape of the vortex plays an underlying role and that the discrete symmetry group of the electronic bands is inherited by the CdGM branches in real space.

### **6.4 Evidence of Multiflux Vortices**

Apart from single-flux vortices with a diameter of ∼ 100 nm, larger vortices could also be observed. Two examples are shown in Fig. 6.10. (a-c) show a vortex that contains two flux quanta. The zero bias conductance map (a) reveals a star shaped STM state pattern in which the arms into the ⟨2¯1¯1⟩ direction are doubled into parallel rays. The arms into the ⟨1¯21⟩ direction remain weak and the conductance in the centre of the vortex is a local minimum (see blue curve in (g)). At higher energy (b), the star arms again split into two arms each totalling four arms in each direction for this vortex. The quasi-classical calculation (a, right panel) confirms that this is the signature of a magnetic vortex with winding number m = 2. Again, the star shaped pattern stems from band II (top right) while band I, with a more isotropic Fermi velocity distribution, would cause a ring (bottom right). Surprisingly, the relative

**Figure 6.10: Multiflux Vortices**: (a-f) dI/dU maps of an m = 2 (a-c) and m = 3 (d-f) vortex. The right panels in (a,d) display the simulated LDOS of band I (bottom) and II (top). (g,h) Single bias spectroscopies (bottom panel) at important locations of the vortex, marked with the corresponding color in the zero bias map (top panel).

magnitude of the LDOS signal from the two bands now seems to be reversed in the calculation. This is in agreement with the experimental observation, where the ring, that is predicted from band I, seems to be overshadowed by the states from band II. An analogous behaviour is found for the vortex with three flux quanta in Fig. 6.10(d-f). Here, the zero bias conductance reveals a star-shaped pattern with three arms in each ⟨2¯1¯1⟩ direction (d) that split into six at higher energies (e). The centre of the vortex now features a local conductance maximum (blue curve in (h)). Both vortices appear larger than the single-flux vortex at eV = ∆<sup>2</sup> (compare Fig. 6.3(h) and Fig. 6.10(c,f)), grow in lateral dimension with the number of confined flux quanta and are not round, but oval-shaped. This might once again be the result of tilted flux lines. The number of confined flux quanta is not limited to small natural numbers: In Fig. 6.11, a giant vortex with m > 10 is shown. It features a large but countable number of flux quanta, because the number of arms that spread from the structure are countable.

Coming back to Volovik's original topological index theorem, one can conclude that for a vortex with anisotropic Fermi velocity, it still holds that the existence/absence of a zero-energy peak in the vortex centre implies an odd/even winding number. However, it is not possible anymore to determine the exact winding number by counting the number of zero-energy crossings in a radial profile. That is because the cylinder symmetry of the problem is removed by the anisotropy in the Fermi

**Figure 6.11: Giant Vortex**: dI/dU maps of a giant vortex with winding number m > 10. The STM states at zero energy (a) extend into the ⟨112⟩ directions with more than 10 arms. At higher energies (b) the arms split like in the vortices with lower m.

velocity. As a consequence, it is not ensured that all diabolic points are crossed in the radial profile, which is crucial for the topological argument. One can, however, now reformulate the topological index theorem for a star-shaped vortex based on the experimental findings and quasi-classical calculations: The number of parallel arms l of a star in the zero-energy LDOS map seems to be identical to the number of diabolic points and thus one can once again use the topological index theorem to state that the winding number of the vortex is given by

$$m = l.\tag{6.3}$$

## **6.5 Vortex Interactions beyond the Ginzburg-Landau Limit**

It is time to address the question, how a stabilization of single-flux and multi-flux vortices in bulk Pb is possible. With a GL parameter of κ ≈ 0.48, Pb is the elemental superconductor closest to the Bogomol'nyi point, κ<sup>c</sup> = 1/ √ 2 ≈ 0.7, but should, according to the standard classification of superconductors, respond to an external magnetic field like a typical type-I superconductor. Apparently, the majority of domains in the intermediate state display the standard behaviour that is characteristic for type-I, i.e. large superconducting and normal regions. During the magnetic field ramp, flux flow is observed. As soon as the magnetic field is kept constant, single vortices could be stabilized, without spontaneous unpinning or drift being observed in the experiment.

This is surprising, because in superconductors of type-I, the ratio of the two characteristic length scales ξ and λ implies an attractive vortex-vortex interaction. That is because the vortex cores overlap before the screening currents see each other and two vortices merge before repulsive interactions come into play. In type-II superconductors, the contrary is the case: Here, the screening currents lead to repulsive interaction between two vortices before the vortex cores overlap. The Bogomol'nyi point with κ = 1/ √ 2 manifests an infinitely degenerate point in the phase diagram where the vortex interaction potential changes sign. In an idealized superconductor, vortices should there behave like non-interacting particles [233–236]. Since the GL theory is only valid near Tc, this standard classification and argumentation is also only valid for temperatures T ≲ Tc. In fact, recent works by Vagov *et al.* on the low temperature extension of the GL theory expect the degeneracy of the Bogomol'nyi point to be lifted below T<sup>c</sup> and give rise to a phase rich transitional region in the (κ, T) plane that widens with decreasing temperature [204]. This is shown in Fig. 6.12. Of special interest is the region of multi-flux vortices right below the transitional domain with single-quantum vortices. The multi-flux vortex region seems to be narrow, but Vagov *et al.* could show in the same publication that a large Fermi velocity anisotropy in two-band superconductors can aid in widening that region further, which opens an alternative path for Pb. To believe that the vortices in this work are in a stable regime due to the temperature dependence of κ only is, however, too naive, because these theoretical calculations neglect the dependence of κ on impurities and the finiteness of sample size.

In order to test whether the above described low temperature effect is sufficient to explain the stabilization of single-flux and multi-flux vortices in bulk Pb, a control experiment at T = 4.3 K was carried out. It turned out, that the stabilization of single vortices within the STM scan frame was substantially harder and only in a single case a stable vortex was found (see Fig. 6.13). This experimental observation suggests that the low temperature plays an important role in the frequent observation of single vortices, but is most likely not solely responsible. As mentioned earlier, in a

**Figure 6.12: Low Temperature Transitional Region**: Low temperature phase diagram of the critical GL parameter κ(T) including the transitional region - κ ∗ 2 , κ<sup>∗</sup> li between standard types I and II. In the region - κ ∗ s , κ<sup>∗</sup> 1 multiquantum vortices are stable solutions. Reprinted figure with permission from [204]. Copyright (2016) by the American Physical Society.

**Figure 6.13: Vortex at** T ≲ Tc: (a,b) dI/dU maps of a magnetic vortex at T = 4.3 K and B = 0 mT. (c) Bias spectroscopies far away from the vortex (black squares) and in the centre of the vortex (blue triangles) in (a). The red line shows a fit of Dynes form to the differential conductance in the superconducting region with a gap size of ∆ = 1.24 meV.

real single crystal of Pb, the GL parameter is also affected by crystal defects, which can at the same time lead to flux pinning as discussed in Chap. 2.1.2.

Extended defects with a size d ∼ ξ are effective pinning centres for flux. An example for such a defect is a screw dislocation, as shown in Fig. 6.14(a). The strong pinning energy scales linearly with the volume V<sup>c</sup> = ξ <sup>2</sup>L, in which the screw dislocation runs through the vortex line, and the condensation energy of Cooper pairs Econd = µ0Hc(0)2/2 [237]. The ends of screw dislocations are clearly visible in STM topographies, but no trend of flux pinning at these sites on the surface was obvious during the experiment. However, at stacking fault tetrahedrons, extended defects that are unique in fcc metals, vortices were regularly pinned. Due to their intricate structure and electronic footprint, a detailed analysis can be found in Chap. 7. In a collective pinning mechanism, defects with a size d ≪ ξ can also effectively trap flux albeit with a smaller scaling ratio. As an example, the weak collective pinning by point defects is depicted in Fig. 6.14(b). In the potential of randomly distributed point defects, the pinning energy scales with <sup>√</sup> Vc. To some extent, the vortex line can sacrifice elastic energy in order to run through parts of the crystal with higher defect concentrations. The length scale L<sup>c</sup> on which the vortex line remains relatively straight is therefore a measure for the pinning strength [237].

All in all, it seems unlikely that the isolated vortices in this work are a result of a high κ alone. In that case, one would expect either a lot more vortices (κ ∗ <sup>1</sup> < κ < κ<sup>∗</sup> li) or domains of all kinds of sizes (κ ∗ <sup>2</sup> < κ < κ<sup>∗</sup> 1 ). Instead, one mostly finds isolated vortices coexisting with µm-large domains. The low temperature does play an important role, as the control experiment at 4.3 K demonstrated. It prevents thermally driven flux creep and puts bulk Pb farther away from the standard type-I regime in the intermediate state so that attractive vortex-vortex interactions are smaller. As a result, lower pinning energies are needed to trap flux on a defect. One can see

**Figure 6.14: Vortex Pinning Mechanisms**: (a) Strong pinning by an extended defect (ddef > ξ), like a screw dislocation: The pinning energy is proportional to the volume, in which the screw dislocation runs parallel through the vortex line. (b) Weak pinning by point defects (ddef ≪ ξ): Small defects can collectively pin flux. The elastic vortex line bends to run through the crystal volume with the most defects. The length scale on which the vortex line remains straight Lc is a measure for the strength of this pinning mechanism. The collective pinning energy only scales like the square root of the total flux line volume. Adapted from Ref. [237].

indications for strong pinning on stacking fault tetrahedrons but not on single screw dislocations. Additionally, impurities can aid the larger defects by a collective pinning mechanism.

### **6.6 Summary**

Although bulk lead (Pb) is classified as a prototypical type-I superconductor, it was demonstrated that the intermediate state hosts single and multi-flux quanta vortices at mK temperatures. As stability criteria for the existence of the vortices, the influence of temperature and defect density on the effective GL parameter κ were discussed. By probing the quasiparticle LDOS inside the vortices with STM and comparing their signature with simulations, it was possible to identify the CdGM states of the two superconducting bands of Pb. The spatial LDOS pattern created by the CdGM states provides a robust method to determine the vorticity or winding number. The robustness is rooted in the fact that the number of diabolic points in the CdGM state spectrum are a topological invariant of the vortex. Simulations of the single and multiflux quanta vortices in the quasi-classical Eilenberger theory, in conjunction with DFT calculated band structure of Pb, support the experimental findings and interpretations. Tilted vortices were useful to demonstrate a low interband coupling for the two superconducting bands, as their respective CdGM state patterns were independently moveable. Additionally, they allowed for a definite assignment of the two gaps ∆ = (1.26 ± 0.02) meV and ∆ = (1.40 ± 0.02) meV to the two superconducting bands of Pb.

# **7 Quasiparticle and FluxTrapping in TopologicalDefects of aTwo-band Superconductor**

Despite being topological defects themselves, vortices in dimension d = 3 are not expected to host topologically non-trivial bound states inside their core [238]. The situation changes when topological crystal defects, i.e. extended defects with dimension D ≥ 1 like dislocations (D = 1) or twin boundaries (D = 2), are added to the mix [239]. In this chapter, the focus lies on stacking fault tetrahedra (SFTs) in Pb(111). These extended defects, that are unique to fcc metals, are efficient pinning centres for flux in the intermediate state, as shown in the differential conductance map in Fig. 7.1. After elucidating their formation and structural properties, their impact on the low energy electronic spectrum is highlighted. Lastly, the interaction of these crystal defects with the CdGM states of a pinned vortex are discussed. This chapter is the logical continuation of Chap. 6 toward the combination of extended

**Figure 7.1: Flux Traps**: Differential conductance map on intermediate state Pb(111) at eU = ∆<sup>2</sup> = 1.4 mV. Areas of low conductance (purple/blue) mark normal conducting regions and areas of high conductance (green/yellow/red) mark superconducting regions. The triangles with a surrounding hexagon of higher conductance are the electronic footprints of SFTs. In two out of three visible SFTs magnetic flux is trapped.

crystal defects and (multiflux) vortices and therefore builds up on findings from the vortex study in Pb.

## **7.1 Stacking Fault Tetrahedra**

In face-centred cubic crystals, the preferred glide plane for dislocations is the (111) plane and the intersections of all available {111} planes form a regular tetrahedron. Therefore, the movement of dislocations is most conveniently explained at the Thompson tetrahedron (TT). As shown in Fig. 7.2, the Burgers vector of a (partial) dislocation can be easily drawn in the TT by the connection of corner points (Latin letters) or face centres (Greek letters).

Perfect dislocations (edge or screw type) are found at the sides of the TT, so the notation for its Burgers vector is the Latin letter for the beginning and end corner point, e.g. BA, which corresponds to the bulk direction <sup>1</sup> 2 [1¯10] in Miller index notation. The unfolded TT in Fig. 7.2(b) serves as a dictionary to translate between Miller index and TT notation. Under applied stress, the glide of a perfect dislocation along a {111} plane is hindered by the hexagonal atomic configuration in the plane and it is energetically favorable for the atoms to move in zig-zag motion via the hcp

**Figure 7.2: Thompson Tetrahedron**: (a) The {111} glide planes in an fcc crystal form a regular tetrahedron - the Thompson tetrahedron (TT). The Burgers vector of (partial) dislocation lines are all found as connections between corner points or face centres of the tetrahedron. (b) The unfolded TT with marked bulk crystal directions serves as a dictionary to translate between Miller indices and TT notation, e.g. BA → <sup>1</sup> 2 [1¯10].

**Figure 7.3: Silcox-Hirsch Mechanism**: (a) The SFT of a frustrated configuration of stair-rod dislocations is built from a triangular Franck partial loop following the dislocation reactions (7.2) and (7.3). Adapted from [240, p. 99]. (b,c) Molecular dynamics (MD) simulations demonstrate that the Silcox-Hirsch mechanism can create an SFT from a triangular (b) or scalene hexagonal platelet of vacancies (c). Reprinted from [241], Copyright (2007), with permission from Elsevier.

hollow site on the (111) plane. Formally, the perfect dislocation dissociates into a pair of *Shockley partial* dislocations:

$$BA \to B\gamma + \gamma A. \tag{7.1}$$

If the stacking fault energy is low enough, the dislocation can travel through the crystal as an extended dislocation, in which a leading Shockley partial and a trailing Shockley partial enclose an intrinsic stacking fault. This means that between the two Shockley partials the stacking sequence in [111] direction is not ABCAB.. but ABABC...

The stacking sequence inside an fcc crystal can also be faulted by the removal or insertion of one close-packed {111} layer of atoms with finite lateral dimension (island). The boundary line of the stacking fault forms a closed loop and is called a *Franck partial* dislocation. The difference to the Shockley partial is that its Burgers vector points perpendicular to the {111} plane, which is faulted, e.g. δD ≡ 1 3 [¯1¯1¯1]. The Shockley partial and the perfect dislocation, whose Burgers vector lie within a {111} plane, are so-called *glissile* dislocations, because they can freely move in the glide plane under applied stress. The Franck partial on the other hand is a *sessile* dislocation and can only move by climb [240, p. 93].

The SFT is a nanocrystal in the shape of a regular or truncated tetrahedron with stacking faults on each of its {111} facets and commonly observed in metals of low stacking-fault energy, like Cu, Ag or Au [242–244]. It is believed to form in the *Silcox-Hirsch* mechanism [240, p. 98-99]: First, vacancies or impurities condense to a triangular platelet, say in the ABC plane, then the bounding Franck partial dissociates on three adjacent {111} facets into a so-called *stair-rod* partial (two Greek letters) and a Shockley partial dislocation (one Greek, one Latin letter):

$$
\delta D \to \delta \alpha + \alpha D \quad \text{(on BCD)}
$$

$$
\delta D \to \delta \beta + \beta D \quad \text{(on ACD)}
$$

$$
\delta D \to \delta \gamma + \gamma D \quad \text{(on ABD)}.\tag{7.2}
$$

The *stair-rod* partials remain in a locked configuration along the triangle edges (on ABC) and repel the glissile Shockley partials due to their parallel Burgers vectors. The Shockley partials then move up the tetrahedron faces and pull up additional stair-rod dislocations on the intersection edges (AD, BD, CD), as shown in Fig. 7.3(a), by the dislocation reactions

$$
\alpha D + D\beta \to \alpha\beta \quad \text{(along CD)}
$$

$$
\beta D + D\gamma \to \beta\gamma \quad \text{(along AD)}
$$

$$
\gamma D + D\alpha \to \gamma\alpha \quad \text{(along BD)}.\tag{7.3}
$$

If the stair-rod dislocations culminate in the tip of the TT (D), then one ends up with a regular SFT, if they stop their movement before, one ends up with a truncated SFT. Atomistic simulations could show that the Silcox-Hirsch mechanism can, under certain circumstances, also develop SFTs from a scalene hexagonal platelet of vacancies [241], which explains the shape of SFTs, observed on Au(111) [244]. The results of these simulations for a triangle and scalene hexagon as the starting vacancy plate are also shown in Fig. 7.3(b,c).

### **7.1.1 Stacking Fault Tetrahedra in Pb(111)**

In the experiments from Chap. 6 and 7, the extended defects which are able to trap flux are interpreted as SFTs based on their size, alignment with the glide planes and electronic footprint. They are most likely created during the hot sputtering process although the dependence of abundance and shape of them on temperature, ion flux and ion energy has not been studied in this work. Figure 7.4(a-d) shows the topographies of four different exemplary SFTs in Pb(111). The white dotted lines correspond to contours of the LDOS as seen in the respective low bias differential conductance maps (e-h) in the superconducting state of the crystal. In (a) and (b) the SFT lies at the surface, which is visible as a scalene hexagonal depression in the topographic image. The apparent negative height is a fraction of one single (111) lattice plane distance. The white contour lines in the [2¯1¯1] and equivalent crystal directions on the surface are in fact the directions in which one would expect projections of the stair-rod dislocations αβ, γα and βγ (cf. Fig. 7.3(a)). The difference in length of these contour lines between (a) and (b) hints towards different evolutionary stages, so different truncation levels, of the SFTs.

By that logic, the SFT in (a) has developed less of a tetrahedron before being truncated by the surface, while the SFT in (b) has developed further, with longer stair-rods. Of course, this is a highly hypothetical interpretation, because the bulk structure is only inferred from a small amount of structural information and the electronic DOS at the surface. Given the fact, that the surface energies in Pb favour vacancy or impurity islands of scalene hexagonal shape in the {111} planes [245], a formation of SFTs like in Fig. 7.3(c) seems most probable. In the topographic images in Fig. 7.4(c,d) a depression like in (a,b) is missing. The SFTs seem to be buried beneath the surface, yet close enough for the LDOS at the surface to confirm their existence. In (c), the SFT seems to bind a non-threading screw dislocation of which both ends are visible at the surface. In (d), no influence on the surface structure by the sub-surface SFT is visible. Unlike the surfacing SFTs, these two SFTs have a rather similar area ratio between the inner scalene hexagon (almost triangle) and the outer hexagon.

**Figure 7.4: Appearance of SFTs**: (a-d) Topographic images of areas containing (truncated) SFTs. The white lines outline where the respective map of differential conductance exhibits clear contours of the LDOS. Surfacing SFTs (a,b) are visible by a hexagonal shaped depression of a fraction of a lattice plane distance. Sub-surface SFTs (c,d) are only detectable through their electronic footprint. (e-f) dI/dU-maps at U = 1.4 meV displaying the corresponding LDOS to the topographies in (a-d). In (f) and (h) flux is threaded through the centre of the defect rendering it normal conducting.

It is possible that these truncated SFTs represent the equilibrium shape acquired universally in the bulk, but a definite answer lacks statistical evidence.

In dI/dU maps at eU = ∆<sup>2</sup> (Fig. 7.4(e-h)) (superconducting state) surface and sub-surface SFTs are easily recognizable by a hexagonal region of ≈ 30 − 150 nm (depending on SFT size) of increased conductance with a cutout, that is usually a scalene hexagon (or almost a triangle), of lower conductance. Rod-like contours of the LDOS can be seen along the sides of the inner triangle/scalene hexagon and along the ⟨1¯21⟩ directions inside the outer hexagon. The former lie within the {111} planes, in which the stacking faults are expected and the latter fall together with the projected directions in which the stair-rod dislocations are expected. In Fig. 7.4 (f,h), flux is threaded through the SFT rendering its inner structure normal conducting (dark purple). The general contour lines of the LDOS, however, remain intact. In Fig. 7.4(e), a standing wave pattern in the green hexagonal (almost round) part of the LDOS is visible, which is absent in the normal conducting state (see Fig. E.1 in Appendix E). In Fig. 7.4(g), this standing wave pattern is either absent or could not be resolved. Here, however, a very clear quasiparticle interference (QPI) pattern is obvious in the inner triangle. Similar particle-in-the-box states have also been found by means of STM in the SFTs of Au(111) and Ag(111) [243, 244] and nanocavities in Pb(111) [245, 246]. The energies at which these states are observed are slightly larger in the SFTs of Au and Ag (E ≈ 1.6 − 4.2 meV) and much larger for the nanocavities (E ≈ 400 meV), which is reasonable considering the size of the respective confinement box. Since stacking-faults are known to have tremendous effect on the electronic band structure, the following section is dedicated to the origin of the QPI pattern and the LDOS contrast around an SFT.

## **7.2 Band Filter Effect and Quantum Well States**

Before turning to the QPI patterns, it makes sense to take a closer look at why the LDOS at the gap edge ∆<sup>2</sup> seems to exhibit such strong local contrast with sharp transitions on the length scale of just a few atoms. As it turns out, this is not limited to the energy ∆2. At eU = ∆<sup>1</sup> the outer hexagon appears lower in differential conductance than the cut-out or the superconducting surrounding but still with the same sharp contour, as Fig. 7.5(a) shows for the SFT from Fig. 7.4(c,g). The source of this contrast reversal in the outer hexagon becomes clear when one looks at the single dI/dU spectra in Fig. 7.5(b) along the line indicated in (a): Inside the inner triangle (points 0-6), the coherence peak at ∆<sup>2</sup> appears to be missing, inside the outer hexagon (points 8-13), the tunnelling amplitude into the coherence peak of ∆<sup>2</sup> exceeds the one for ∆<sup>1</sup> and outside the hexagon (points 15-25), the spectra regain the characteristic coherence peak heights for the Pb(111) surface. This behaviour is best seen in a plot of the fitted amplitude σ<sup>i</sup> of the BCS curve for the gap ∆<sup>i</sup> as a function of spectroscopic location, like in Fig. 7.5(c). In the bottom panel, the ratio A = σ<sup>2</sup> σ<sup>1</sup> and the asymmetry χ = σ1−σ<sup>2</sup> σ1+σ<sup>2</sup> are plotted. While the tunnelling amplitude into the two coherence peaks of Pb was demonstrated to be dependent on the tunnelling direction with respect to the crystal orientation [39], here, only vertical tunnelling into (111) oriented planes is expected above the SFT<sup>1</sup> . Together with the fact that the superconducting properties cannot change on a length scale much

**Figure 7.5: Band Filter Effect**: (a) dI/dU map (eU = ∆1) of SFT from Fig.7.4(c,d). Black crosses mark single spectra location and the red arrow indicates the direction in which the points were sequentially probed. (b) Two-gap BCS fits to the spectra obtained along the line marked in (a). The single spectra are offset in y-direction for visibility. (c) Top panel: Scaling factor (amplitude) σ1/σ<sup>2</sup> of the BCS curves with ∆1/∆<sup>2</sup> along the line marked in (a). Bottom panel: Amplitude ratio A = σ2/σ<sup>1</sup> and asymmetry χ = (σ<sup>1</sup> − σ2)/(σ<sup>1</sup> + σ2) along the same line.

<sup>1</sup> Even in the presence of screw dislocations, the distorted planes around the dislocation retain the character of (111) planes because there is no stacking fault.

shorter than the coherence length, it is most likely that the seen effect is caused by a constraint imposed on the movement of normal conducting electrons in the system. It seems that the stacking faults of the SFT act as a band filtering energy barrier that traps QPs of band I with the energy ∆<sup>1</sup> inside the tetrahedron and hinders QPs of band II with energy ∆<sup>2</sup> from entering it. As a consequence, a QPI pattern of QPs from band I<sup>2</sup> is formed inside the tetrahedron (see Fig. 7.5(a)) and standing waves of QPs from band II are formed outside the tetrahedron (see Fig. 7.4(e)) as a result of their backscattering. What the authors of Ref. [247] claim to be "enhanced superconductivity" around nanocavities in Pb(111) could also be explained by the band selective scattering of QPs with energy ∆<sup>1</sup> and ∆2, which would not be resolved in their experiment due to insufficient energy resolution. What is puzzling, is that the normal conducting state, so far, showed no signs of this band filter effect.

**Figure 7.6: Scattering Vectors**: (a) dI/dU map (eU = ∆2) of the SFT from Fig. 7.4(a) reveals standing wave fronts in the green region oriented toward the ⟨2¯1¯1⟩ and the ⟨1¯21⟩ directions. (b) FFT of the region marked in (a). The scattering vector (red arrow) has the length |q∥| = 0.78 π/a in the [2¯1¯1] bulk ([11] surface) direction. (c) Corresponding nesting vector q<sup>∥</sup> = (0.65, 0.375) π/a on the Fermi surface of band II responsible for the scattering vector in (b). (d) dI/dU map (eU = ∆1) of the SFT from Fig. 7.4(c) reveals QPI pattern inside the inner triangle. (e) Analogous to (b): the scattering vector is |q∥| = 0.26 π/a in the [11¯2] direction. (f) Analogous to (c): The corresponding nesting vector on the FS of band I has q<sup>⊥</sup> ≈ 0.333 π/a and q<sup>∥</sup> = (0, 0.26) π/a. In (c,f) k<sup>x</sup> and k<sup>y</sup> point in the [1¯10] and [11¯2] direction. The axes in (b,e) are slightly misaligned to these directions (see real space directions in Fig. 7.4(a)).

<sup>2</sup> The QPI pattern seen in the same defect at eU = ∆<sup>2</sup> in Fig. 7.4(g) is caused by the oscillation of the coherence peak amplitude at ∆<sup>1</sup> and the sensitivity of the measurement to its high energy flank at ∆2.

In a local STS measurement, one gains exceptional real space resolution for the sacrifice of momentum information, because the pure eigenstates in real space are the Fourier transformed eigenstates in k-space. Techniques like angle-resolved photoemission spectroscopy (ARPES) in turn probe the pure momentum eigenstates that contain no local information. In the presence of crystal impurities, however, these two pictures are mixed because the translation invariance is broken. Especially in inhomogeneous superconductors like Bi2212 [248, 249], researchers developed the technique of Fourier-transform STM, where the probed LDOS of the QPs could be linked to scattering vectors on the Fermi surface (FS) and ultimately explain discrepancies between STM and ARPES studies. The coherent scattering on a defect leads to oscillations of the LDOS in real space with the period λ = 2π/|q|, where q = k − k ′ is the scattering vector.

The Fast-Fourier transform<sup>3</sup> of the QPI pattern marked in Fig. 7.6(a) in the outer hexagon reveals the scattering vector with length |q∥| = 0.78 π/a in the [2¯1¯1] direction (see (b)). This vector connects opposite sides of the Fermi surface of band II (along U−Γ−U andK−Γ−K) where the DOS is relatively small but this is overcompensated by the connection of very large portions of the FS by this single vector, as shown in (c). The same technique can be applied to the QPI pattern in Fig. 7.6(d) at energy ∆1. One finds the scattering vector with |q∥| = 0.26 π/a in all ⟨2¯1¯1⟩ and ⟨1¯21⟩ directions in the (111) plane, as shown in (e). The corresponding 3D scattering vector with the same q∥, shown in (f), connects opposite sides of the smaller necks in the FS of band I (along U − X − U) in the Γ − X direction. Also here, the connected FS parts are almost parallel to each other and have a relatively high DOS. Because the QPI pattern in (d) is limited to a closed box (the inside of the SFT), one needs to be a bit more careful in the interpretation. That is because the boundaries of the box impose another (larger) periodicity on the electronic wavefunctions on top of the lattice periodic crystal potential.

The Schrödinger equation of a free particle in an equilateral triangle (2D) box has an analytical solution [250, 251] that has been experimentally demonstrated by the authors of Ref. [243] in the STM study of SFTs on Ag(111). Ag(111) has a surface state of quadratic dispersion and SFTs with a triangular ground plane making it an optimal system to study the QWSs of the triangular 2D box. The electrons in this surface state are a quasi-free 2D electron gas with constant DOS so that the interference pattern inside the SFT is solely determined by the box geometry and energy. On Pb(111), such a surface state is absent, because band II does not have a bulk band gap in this direction. In nanocavities of Pb(111) [245, 246], the QWSs are still dictated by the box geometry because the box size is of the same order as the

<sup>3</sup> A Hanning window function was used in the Fast-Fourier transforms.

**Figure 7.7: Free particle in an Equilateral Triangle Box**: Analytical solution for the LDOS of a free electron with quantum numbers (p, q) = 2 3 , 17 <sup>2</sup> 3 in an equilateral triangle following Ref. [250].

lattice constant and therefore overshadows the bandcharacter of the electrons. The QPI in an SFT of superconducting Pb(111) like in Fig. 7.6(a) should in principle be viewed as the result of Bogolyubov QPs just above the gap energy that inherit the Bloch character of electrons<sup>4</sup> at the Fermi surface and scatter at the stacking faults of the tetrahedron. Therefore, they should scatter predominantly with scattering vectors that connect parallel parts of the FS with high DOS and indeed this is what was seen in Fig.7.6.(d-f). Still, the interference pattern shows striking resemblance to the QWS obtained in the model of a free electron gas inside a triangular 2D box, that is displayed in Fig. 7.7. Considering the fact that the QWSs are not (or not clearly) visible in the defect from Fig. 7.6(a), it seems that in Fig. 7.6(d) the SFT/triangular box has just the right size to support QWSs with a periodicity in real space that corresponds to strong scattering vectors.

In order to shed more light on the nature of the interference pattern, it should first be made sure that the same pattern is seen at opposite bias voltage which is to be expected for the near-gap QPs of a superconductor. Secondly, the evolution of the interference pattern with energy could be measured in order to learn more about the dispersion of the QPs. And lastly, a larger number of SFTs should be investigated to make concluding remarks on the influence of the box size and shape.

<sup>4</sup> In this case, the box size is much larger than the periodicity of the lattice potential making a description in the basis of Bloch states reasonable.

## **7.3 Electronic Edge State**

Stacking faults have already proved to be efficient strong pinning centres for flux, at least in high-temperature superconductors with columnar defects [252, 253]. Because a tendency toward flux pinning at SFTs was seen in this work, it was possible to pin a single-flux vortex inside an SFT.

In Fig. 7.8(a), the zero-bias conductance map around the defect from Fig. 7.4(c) features the star-shaped CdGM states characteristic for a single-flux vortex. The central line of the vortex (centre of CdGM state pattern) lies inside the SFT as a comparison to the dI/dU map at eU = ∆<sup>2</sup> (e) clarifies. The star arms in the ⟨2¯1¯1⟩ directions, which stem from band II, are shorter than before. The states from band I, so the weak arms in ⟨1¯21⟩ directions and the higher energy CdGM states (b-c), that produce a ring in the real space pattern, do not seem to cross the triangular outline of the SFT. A logical explanation would be that this effect is also caused by the band filtering of the stacking faults. Although not well resolved, in (b,c) the high-energy CdGM states seem to produce a QPI pattern similar to the one observed for the gap edge QPs. One additional feature is striking in the map in (a): In the presence of the pinned single-flux vortex, an additional state, which is localized on the edge of the

**Figure 7.8: Single-flux Vortex in SFT**: dI/dU maps of the SFT from Fig. 7.4(c) at B = 6 mT: a single-flux vortex is threading the tetrahedron. (a) At U = 0.0 mV the CdGM state pattern of the vortex is less extended (shorter star arms) and along the SFT edge a zero-energy edge state emerged. (b-e) At U = 0.4, 0.8, 1.2, 1.4 mV, the regular behaviour of CdGM states in Pb is seen. The CdGM states also seem to produce a QPI pattern inside the tetrahedron.

**Figure 7.9: Edge State**: (a) dI/dU map of same area as in Fig. 7.8 after decreasing the magnetic field to B = −2 mT. The vortex moved further up but remains inside the SFT and the edge state prevails. (b) Line profile of the differential conductance as function of bias voltage along the white line marked in (a). The CdGM states do not cross zero energy. At the position of the SFT edge a very localized peak centred around zero energy emerges.

SFT, appeared. This state is absent in the purely superconducting phase (see Fig. E.2 in Appendix E).

Even after slightly shifting the vortex line by a reduction of the applied magnetic field from B = 6 mT to −2 mT, the edge state persists, as demonstrated in Fig. 7.9(a). A closer look at the differential conductance spectrum along the white dotted line (b) makes it clear that the edge state is distinct from the CdGM branch in the spectrum. The latter does not move considerably in energy along the line profile, the former appears exactly at the triangular contour of the SFT where the stacking fault line is assumed to be.

In order to continue the discussion of inherent differences between the CdGM states and the edge state on the basis of physical quantities, the localisation and "lifetimes"<sup>5</sup> or spectral bandwidths of these states were compared. This is shown in Fig. 7.10.

<sup>5</sup> The spectral bandwidth of CdGM states around zero energy will for reasons of readability be associated with a "lifetime", although, this interpretation is strictly not allowed. In a measurement at a single tip position, several low energy CdGM states are expected to overlap. STS can never distinguish between single life-time broadened states and bands if the energetic spacing of states is below the resolution capability of the instrument. The general conditions for each single measurement are, however, the same, so the deduced spectral widths will still only differ if the properties of the states themselves are different.

**Figure 7.10: Localisation and Lifetimes**: (a-c) Zero bias dI/dU maps of normal single flux vortices trapped on an SFT (a,b) (vortex from Fig. 7.8, 7.9) and of an anomalous vortex (b) (vortex from Fig. 6.5) with marked cross sections C1-18 (lines) and point spectroscopies P1-3 (circles). The localisation of the zero energy edge state (d), the zero energy CdGM state of band II (f) and the zero energy CdGM state of band I (h) was determined through a fit of the average cross section dI/dU(U = 0, x) to a Gaussian. The corresponding lifetimes were determined by a fit of the point bias spectroscopies to a Lorentzian (e,g,i).

In the normal vortex pinned on an SFT (a) the localisation of the edge state can be extracted from cross sections of the dI/dU signal at U = 0 mV across the edge. To enhance the signal-to-noise ratio, several cross-sections (C1-7) were averaged. This averaged line profile was fitted to a generalized Gaussian of the form

$$f(x) = \frac{C}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(x-x\_0)^2}{2\sigma^2}\right) + b + cx \tag{7.4}$$

with the characteristic width σ. The lifetime τ = 1/Γ was determined from a single point spectroscopy on the edge (P3) shown in the pinned vortex in (b). The dI/dU(U) curve was reduced by a constant background and then the part around U ≈ 0 was fitted to a Lorentzian of the form

$$g(U) = \int\_{-\infty}^{\infty} \mathrm{d}U' \frac{C}{(U'-U\_0)^2 + (\Gamma/2)^2} \frac{\partial n\_F(U'-U)}{\partial U} \tag{7.5}$$

with the characteristic scattering rate Γ = 1/τ . For the CdGM states, the fact that the states from the different bands can be probed independently, especially in an anomalous vortex (c), was used. The localisation and lifetime of the zero energy CdGM states of band II (star pattern) were extracted from the cross-sections C8-17 (f) and the point spectroscopy P1 (g). The results for the zero energy CdGM states of band I (point at P2) are shown in (h) and (i). The respective localisation widths σ and scattering rates Γ from the fits are indicated in the figure. Even from this small set of data, significant differences between edge state and the CdGM states

vortex showing no sign of an edge state at the triangular contour of the SFT (white dotted line). (d-f) Surfacing SFT with a pinned single-flux vortex showing a clear edge state at zero bias. (a,d) Topographic images of the area containing the SFT. (b,e) dI/dU maps at U = 1.4 mV of the flux-penetrated SFT. (c,f) Corresponding dI/dU maps at U = 0 mV.

**Figure 7.12:** m = 2 **Vortex in SFT**: dI/dU map at U = 0.0 mV of an m = 2 vortex trapped inside the sub-surface SFT from Fig. 7.4(d). No clear edge state is visible.

become apparent. The edge state has a much smaller lifetime than the zero energy CdGM states of band I, it is more comparable to the zero energy CdGM states of band II. It is, however, much more localized than the CdGM states of both bands as a comparison of the deduced localisation widths σ shows. From this analysis, a clear one-to-one correspondence between the edge state and the CdGM states is not seen. The fact that a sudden change from normal conducting region to superconducting region is seen at the edge of this SFT is intriguing in the sense that this edge state may play an important role for the flux screening.

Comparing two SFTs with a pinned single-flux vortex, one buried beneath the surface (Fig. 7.11(a-c)) and one at the surface (Fig. 7.11(d-f)), indicates that the edge state (marked with white arrows) is only visible for surfacing SFTs. This would not be too surprising, considering the fact that the edge state is also strongly localized in the lateral direction. Hence, it might be present but not detectable on the surface for buried SFTs. Of course, for a convincing argument a higher number of SFTs with pinned vortices would have to be analysed.

The pinning of an m = 2 vortex inside the defect from Fig. 7.4(d) has also been achieved. In the zero-bias conductance map in Fig. 7.12, the characteristic footprint of the m = 2 vortex is seen, yet no edge state is visible, where the SFT edge is seen in the map at eU = ∆<sup>2</sup> (cf. Fig. 7.4(h)). Whether the invisibility of the edge state in this case is due to the even vorticity of the vortex, the fact that this defect is buried deeper beneath the surface or has an entirely different reason, cannot be concluded from the available data.

While a topological origin of this state is far from certain, none of the here reported properties are in conflict with it. In fact, some of them constitute necessary conditions for a topological edge state. In summary, the following properties of the SFT edge state were found:


## **7.4 Summary**

By hot sputtering a Pb(111) single crystal with keV argon ions and subsequent annealing, SFTs were formed. In the superconducting state, these extended crystal defects have tremendous influence on the LDOS of QPs. They are responsible for a memorable LDOS pattern at the gap energy that cannot be explained by structural deformations alone. A large change in the electronic band structure is expected at the stacking faults that surround the tetrahedral nanocrystal. This boundary seems to act as a band selective QP filter. It traps QPs of one band inside the SFT in the form of quantum well states (QWSs) and prevents the QPs of the other band from entering. Since analogous measurements in the normal state do not show these effects, the SFT may potentially constitute a Cooper pair box. Upon pinning vortices on SFTs, the transition between normal conducting region and superconducting region is in certain cases reduced to a length scale much lower than the known bulk coherence length. Concurrently with this, a zero-energy edge state emerges at the boundary of the SFT. A series of systematic tests with the purpose of fathoming the origin of this state could not reject a topological cause.

# **8 Conclusion andOutlook**

In the scope of this thesis, the full potential of a mK scanning tunnelling microscope (STM) was used to study three very different types of superconducting materials that are constant subject of fundamental research: A granular superconductor, a conventional multi-band superconductor and unconventional cuprate superconductors. In the presence of localized spins, vortices or crystal defects the elastic part of the tunnelling current revealed a large variety of interesting quasiparticle (QP) effects that include Yu-Shiba-Rusinov (YSR) states, Caroli-de Gennes-Matricon (CdGM) states and interband couplings. Outside of the superconducting gap, inelastic tunnelling events contributed significantly to the tunnelling current. A special focus was given to antiferromagnetic spin fluctuations and potential plasmon or Higgs modes near a superconductor-to-insulator transition (SIT). In the analysis part of the thesis, the deeper understanding of CdGM states and observed inelastic excitations required quasiclassical calculations in the framework of 3D Eilenberger equations and an advanced deconvolution algorithm that were independently programmed by the author.

In quest of the bosonic glue in cuprate superconductors, we demonstrated how the effective Eliashberg function can be extracted from IETS data in the superconducting state. The extracted boson spectral function of these d-wave superconductors is in good agreement with the expected integrated spin spectrum, which validates the spin-fermion model. As momentum information is usually lost in a tunnelling experiment, an essential ingredient in the success of this method was the dominant influence of hot-spots on the Fermi surface on the superconducting density of states. The real space imaging power of the STM could be used further to examine the position dependence of the bosonic excitation and potentially even perform QPI measurements at the right energies. The minimum energy of the inelastic scattering channel is expected to vary locally due to the inhomogeneity of ∆.

We put together the first STM study of granular aluminium (grAl) (that we know of) and published the results in Physical Review B [146]. The systematic study of samples with different oxygen content could confirm many microscopic properties of these superconducting films that were previously derived from bulk measurements, like electronic transport or optical spectroscopy. These include a homogeneous superconducting gap, as well as a gap enhancement and grain decoupling that scale with the oxygen concentration. In addition, on individual grains in the highresistivity film, the tunnelling electrons could easily excite a long-lived bosonic mode. It is likely that this mode corresponds to phase (or amplitude) fluctuations of the superconducting gap near the SIT. Perhaps the most important result in the high-resistivity grAl film was the observation of YSR states. The careful analysis of position dependent tunnelling spectra hints towards localized spins that reside in the oxide barriers between the aluminium grains. If this is indeed the case, then this has important implications for many Al/AlO<sup>x</sup> based superconducting devices as spin scattering can potentially limit their performance. In principle, these localized spins should also appear in the natural oxide of aluminium. First attempts to find YSR states in the natural oxide of a polycrystalline aluminium film could, however, not confirm them. The measurements should be repeated on a thicker oxide before final conclusions can be drawn. As the tunnelling barrier increases with a thicker surface oxide layer, different experimental techniques may be required.

It is also worthwhile to revisit supposedly well understood conventional superconductors like lead (Pb) with state-of-the-art STM setups as ultra-low temperatures and better energy resolution can provide new insights into the microscopic understanding of the superconducting properties. In this work, it was demonstrated that the intermediate state of a bulk Pb single crystal can naturally host stable vortices with single and multiple flux quanta at mK temperatures. We thus have to rethink our traditional classification of superconductors into type-I and type-II on the microscopic scale and at temperatures far below the critical temperature. Using the exceptional energy resolution of a 25 mK STM, it was possible to study the interplay of CdGM states from two different bands inside a single vortex line. Additionally, the long theoretically predicted LDOS patterns of CdGM states in multiflux vortices could for the first time be confirmed experimentally and their immediate connection to the vortex winding number was demonstrated. While the CdGM states of a vortex in an s-wave superconductor are not considered to be topologically protected bound states, topologically non-trivial boundary states can in principle exist at topological defects. A combination of the multiflux vortices with extended crystal defects seemed like the obvious next step.

We saw that stacking fault tetrahedra (SFTs) have a detrimental effect on the local quasiparticle DOS of Pb. Apparently, the weak interband coupling of the two superconducting bands leads to band selective scattering processes at the SFT boundary and confinement effects. Stacking faults are known to have a relatively large effect on the electronic band structure of the normal state, yet it remains an open question which role superconductivity plays in this case. The SFTs can also act as pinning centres for multiflux vortices. As it turns out, the CdGM states are also affected by the SFT boundary and in certain cases, a pinned single-flux vortex lead to a zero-energy edge state localized on the SFT boundary. Whether this state is topological in the quasiparticle sense remains to be seen, but the present findings do not yield a categorical argument against it.

Superconducting phenomena often span very different classes of superconductors. One example is the SIT that is of central importance in various high temperature superconductors (HTSs), but also in granular superconductors. Another example is multiband superconductivity which is widely accepted to be present in ironbased superconductors and the elementary superconductor Pb but has been mostly overlooked up until the discovery of two-band superconductivity in MgB<sup>2</sup> in 2001. This thesis shows that it is worthwhile to revisit superconductors that have been around for several decades and that we can, when equipped with state-of-the-art experimental and theoretical techniques, reveal previously hidden properties. When being open for surprises, we can often find them in places where we did not expect them. Of special noteworthiness are the findings of localized spins in grAl with important implications for the use of Al/AlO<sup>x</sup> based superconducting devices and the existence of stable (multiflux) vortices in bulk Pb, which provide optimal conditions to study vortex core states in combination with multiband effects. Especially with the discovery of an intriguing edge state, that emerged where vortex states were perturbed by extended crystal defects, we hope to initiate further research towards the realization of a Majorana zero mode (MZM).

# **Bibliography**


# **Acronyms**

**AC** Alternating current **ADC** Analog-to-digital converter **AES** Auger electron spectroscopy **AFV** Antiferromagnetic ordering vector **ARPES** Angle-resolved photoemission spectroscopy **BCS** Bardeen-Cooper-Schrieffer **BdG** Bogoliubov-de Gennes **BEC** Bose-Einstein-Condensate **Bi2212** Bi2Sr2CaCu2O8+<sup>δ</sup> **BSCCO** Bismuth strontium calcium copper oxide **BZ** Brillouin zone **CdGM** Caroli-de Gennes-Matricon **D-STM** Dilution scanning tunnelling microscope **DAC** Digital-to-analog converter **DC** Direct current **DFT** Density functional theory **DOS** Density of states **FFT** Fast Fourier transform **FIM** Field ion microscopy **FPGA** Field programmable gate array **FS** Fermi surface **FWHM** Full width at half maximum

**GL** Ginzburg-Landau **grAl** Granular aluminium **HTS** High temperature superconductor **IETS** Inelastic tunnelling spectroscopy **INS** Inelastic neutron scattering **IQMT** Institute for Quantum Materials and Technologies **irrep** Irreducible representation **JJ** Josephson junction **JT-STM** Joule-Thomson scanning tunnelling microscope **LBCO** Lanthanum barium copper oxide **LDOS** Local density of states **LEED** Low energy electron diffraction **MBE** Molecular beam epitaxy **MIT** Metal-to-insulator transition **MO** Magneto-optical **MZM** Majorana zero mode **NEG** Non-evaporable getter **NG** Nambu-Goldstone **NIN** Normal metal-insulator-normal metal **NIS** Normal metal-insulator-superconductor **NMR** Nuclear magnetic resonance **QMS** Quadrupole mass spectrometer **QP** Quasiparticle **QPI** Quasiparticle interference **QWS** Quantum well state **RC** Real-time controller **REBCO** Rare earth barium copper oxide

**RMS** Root mean square

**SCBA** Self-consistent Born approximation

**SEM** Scanning electron microscopy

**SFT** Stacking fault tetrahedron

**SIS** Superconductor-insulator-superconductor

**SIT** Superconductor-to-insulator transition

**SQUID** Superconducting quantum interference device

**STM** Scanning tunnelling microscope

**STS** Scanning tunnelling spectroscopy

**TMP** Turbomolecular pump

**TSP** Titanium sublimation pump

**TT** Thompson tetrahedron

**UHV** Ultra-high vacuum

**WKB** Wentzel–Kramers–Brillouin

**Y123** YBa2Cu3O6+x

**YBCO** Yttrium barium copper oxide

**YSR** Yu-Shiba-Rusinov

**ZFC** Zero field cooling

# **Notation**

Throughout this thesis, the SI unit system and the following notations were used:


#### **List of physical constants**


# **List of Figures**




# **List ofTables**


# **Appendices**

# **A Deconvolution using theGold Algorithm**

On a very general basis, the spectroscopic output signal of a system y with a certain response function h can be written as the discrete convolution

$$y(\omega\_i) = \sum\_{j}^{i} x(\omega\_j) h(\omega\_i - \omega\_j) + n(\omega\_i),\tag{A.1}$$

where x is the input signal, n is arbitrary noise and ω<sup>i</sup> are discrete increments of frequency. In Chap. 4, y ↔ σ i , x ↔ g 2χ ′′ and h ↔ σ el. One can write Eq. (A.1) in matrix form [254]:

$$y = Hx + n.\tag{A.2}$$

A least square solution of these coupled linear equations is found by minimization of the functional [254]

$$\|H\hat{x} - y\|\_2 \, , \tag{A.3}$$

where xˆ = A<sup>−</sup>1Hy is the unconstrained least squares estimate of the input signal and A = H<sup>T</sup> H is called the Toeplitz matrix. This solution is not definite, because the problem of finding x with knowledge of H and y is an ill-posed problem. There exist many solutions to this system of equations and small errors in assumption of H or noise cause large oscillations in the extracted solution. Regularization techniques are necessary to confine the possible solutions to a smaller subset. Common regularization techniques concentrate on modelling the noise. In the problem from Chap. 4, the shape of the noise is unknown but even worse, the input signal x is only a vague guess. Therefore, larger errors due to the uncertainty of H are expected to influence the deconvolution to find x. A physical constraint that drastically reduces the amount of possible solutions is the required positiveness of y, x and therefore also the positive definiteness of the Toeplitz matrix A.

The Gold algorithm [141, 254] is one of the few deconvolution algorithms that makes use of this constraint and always yields positive solutions if the input is positive. In the classic Gold algorithm, the solution is refined iteratively by the procedure [254]

$$x^{(n+1)}(\\\omega\_i) = \frac{y'(\omega\_i)}{\sum\_{j=0}^{M-1} A\_{ij} x^{(n)}(\omega\_j)} x^{(n)}(\omega\_i). \tag{A.4}$$

This solution converges to the least square estimate in the constrained subspace of positive solutions and is frequently used in the analysis of γ-spectroscopy data where detector counts are necessarily positive. The algorithm also decomposes the solution into a series of δ-peaks which is not ideal for the problem at hand, but can pronounce the seeked bosonic mode. Since the dimension N of A takes at least the number of points in the spectrum, it is large and the product of Toeplitz matrix and vector x in Eq. (A.4) is done by the FFT. In fact, the kernel function is far from periodic and does not naturally go to zero for ω → ∞ which is why it has to reasonably padded, meaning it is first extrapolated and then smoothly brought to zero. The total deconvolution range then spans a range much larger than the range where the kernel is unequal to zero. A typical total size of the Toeplitz matrix is N = 12000 due to the extensive padding that is required. The FFT thus reduces the computation time for one step by a factor of

$$\frac{N^2}{\frac{1}{2}N\log\_2 N} \approx 1770.\tag{A.5}$$

# **B SupplementaryGrAlData**

In this appendix, supplementary data supporting the conclusion from Sec. 5.4 is presented. Fig. B.1 shows a spectrum containing out-gap peaks indicating inelastic excitations of the tunnelling electron, just like in Fig. 5.9 of the main text. This spectrum was obtained at a different tip position and slightly milder tunnelling conditions, so smaller overall conductance at the feedback condition (U = 4 mV, I = 70 pA). It was thoroughly checked that periodic noise in the tunnelling contact is not causing the equidistant peaks in the spectrum. For good measure, the bias step size was changed (while keeping the integration time for each bias value equal) compared to the spectrum of Fig. 5.9. The spectrum shown here is an average of 1000 single spectra. For the analysis of this spectrum, the model from Eq. (5.2) of the main text was fitted to the data. Tab. B.1 shows the found parameters, when fitting first the negative bias side for the energies and then imposing them for the positive side. Again, with five excitations, the spectrum is qualitatively well reproduced (see red line in Fig. B.1). Especially the ω<sup>n</sup> for n = 2, 3, 4 are very similar to the ω<sup>n</sup> found for the data presented in the main text. From this analysis, ω<sup>1</sup> is slightly larger than 1.35 ∆.

The deconvolution by the FFT and Gold algorithm, however, reveals peaks in the boson spectral density α <sup>2</sup>F(ω) at multiples of ωexc = 1.35 ∆, as indicated by the grey dashed lines in the bottom panel of Fig. B.2. Solely the fourth excitation is found at higher energies than expected. The trend of increasing amplitudes for harmonics of n < 3 and subsequent decrease for harmonics of n > 3 is in good agreement with what was interpreted into the boson spectral density obtained in the main text.

**Table B.1: Fit Parameters for Bosonic Excitations**: Fitted parameters of Eq. (5.2) used to model the spectrum in red in Fig. B.1. n is the number of the bosonic excitation, ω<sup>n</sup> its energy divided by ∆ = (250 ± 3) µeV and a + <sup>n</sup> and a − n are the intensities for positive and negative bias range.

**Figure B.1: Inelastic Excitations**: Analogon to Fig. 5.9, recorded at a different tip position. Here, the left side was used for the determination of the five excitation energies. Similar energies as in the main text are found (see Tab. B.1).

**Figure B.2: Boson Spectral Density**: Analogon to Fig. 5.10 using the spectrum in Fig. B.1 as a source. Grey dashed lines indicate the energies n · ωexc for n = 1, 2, 3, 4 and ωexc = 1.35 ∆ like in the main text.

# **C 3D Eilenbenberger Equations**

Gert Eilenberger [255] and by a slightly different approach Larkin and Ovchinnikov [256] could show that the Gorkov equations of motion for normal and anomalous Green's functions can be transformed into a set of simpler transport-like equations for quasi-classical Green's function propagators, that hold for k<sup>F</sup> ξ ≫ 1 [257]:

$$-\hbar v\_F \nabla \hat{g}(\mathbf{r}; \mathbf{p}\_F, i\epsilon\_n) = \begin{bmatrix} \left(i\epsilon\_n + v\_F e\mathbf{A}(\mathbf{r}) & -\Delta(\mathbf{r}, \mathbf{p}\_F) \\ \Delta^\dagger(\mathbf{r}, \mathbf{p}\_F) & -i\epsilon\_n - v\_F e\mathbf{A}(\mathbf{r}) \end{bmatrix}, \hat{g}(\mathbf{r}; \mathbf{p}\_F, i\epsilon\_n) \right). \tag{C.1}$$

The quasi-classical Green's function propagator

$$\hat{g} = \begin{pmatrix} g(\mathbf{r}, \mathbf{p}\_F, i\epsilon\_n) & f(\mathbf{r}, \mathbf{p}\_F, i\epsilon\_n) \\ -f^\*(\mathbf{r}, -\mathbf{p}\_F, i\epsilon\_n) & g^\*(\mathbf{r}, -\mathbf{p}\_F, i\epsilon\_n) \end{pmatrix},\tag{C.2}$$

depends on spatial coordinate r, crystal momentum p<sup>F</sup> = ℏk<sup>F</sup> and energy ϵn, and must satisfy the normalization condition gˆ <sup>2</sup> = ˆ1. g and f are normal and anomalous quasi-classical Green's function propagators and ϵ<sup>n</sup> are fermionic Matsubara frequencies. Incorporated in this formalism is the self-consistent calculation of pair potential ∆(r) and A(r), which is essential for the description of a superconductor hosting vortices as these two fields are dependent on each other. Although this approach is numerically much more feasible than direct diagonalization of a BdG Hamiltonian when treating inhomegeneties, it is not necessary to solve the problem completely self-consistently in this work. The reason for that is, that the experimental results from Chap. 6 clearly show, that, even though the sub-gap states in the vortex display highly anisotropic behaviour, the recovery of ∆(r) is isotropic in the plane around the vortex core. The local pair potential is assumed to have s-wave symmetry and has therefore been modeled by

$$\Delta(r, p\_F) = \Delta(r) = \left(\Delta\_0 \tanh\frac{r}{\xi} + \Theta(z)W \tanh\frac{z}{a}\right)\left(\frac{x+iy}{r}\right)^m,\tag{C.3}$$

with Θ being a Heaviside step function, ∆<sup>0</sup> the maximum gap size, ξ the coherence length, a the lattice constant, m the winding number of the vortex and W the work function. The second term ensures, that quasiparticles travelling to the surface are decaying into the vacuum with the right damping factor. Since an isotropic ∆ implies an isotropic in-plane current density, the magnetic field profile around the vortex is described by a vector potential of the form [258]

$$\mathbf{A} = A(r, z)\hat{\mathbf{e}}\_{\varphi} \,, \tag{\text{C.4}}$$

$$A(r,z) = \frac{m\Phi\_0}{2\pi\lambda^2} \int\_0^\infty \mathrm{d}k \frac{J\_1(kr)}{k^2 + \lambda^{-2}} S(k,z) \,, \tag{C.5}$$

$$S(k, z) = \begin{cases} \frac{\kappa}{k + \kappa} \mathbf{e}^{-kz} & z > 0 \\ 1 - \frac{k}{k + \kappa} \mathbf{e}^{\kappa z} & z \le 0 \end{cases} \tag{C.6}$$

that is cylinder symmetric in the bulk and deviates from the bulk value near and above the surface. Here, κ = √ k <sup>2</sup> + λ<sup>−</sup>2, λ is the magnetic penetration depth which was chosen to be <sup>λ</sup> <sup>=</sup> ξ/<sup>√</sup> 2 and J1(x) is a Bessel function of first order.

Nils Schopohl and Kazumi Maki [259] could show that the Eilenberger equations can always be solved along a characteristic line and that the solution is universal for every point along this line. The line simply has to be parallel to the Fermi velocity vector v<sup>F</sup> . Points on this line are then characterized by the variable X and two *impact parameters* Y and Z.

$$
\sigma(X) = X\hat{\mathfrak{u}} + Y\hat{\mathfrak{v}} + Z\hat{\mathfrak{w}},\tag{C.7}
$$

$$
\hat{\mathbf{y}} = x\hat{\mathbf{x}} + y\hat{\mathbf{y}} + z\hat{\mathbf{z}}\tag{C.8}
$$

Transformation from the mobile frame of reference (uˆ, vˆ, wˆ), where uˆ ∥ v<sup>F</sup> , to the fixed coordinate system (xˆ, yˆ, zˆ) is done by chaining rotation matrices using Euler angles:

$$
\begin{pmatrix} v\_x \\ v\_y \\ v\_z \end{pmatrix} = R\_{zyz} \begin{pmatrix} v\_X \\ v\_Y \\ v\_Z \end{pmatrix}, \tag{C.9}
$$

$$\begin{aligned} R\_{zyz} &= R\_z(\eta) R\_y(\chi - \pi/2) R\_z(\psi) \\ &= \begin{pmatrix} \cos \eta & -\sin \eta & 0 \\ \sin \eta & \cos \eta & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \sin \chi & 0 & \cos \chi \\ 0 & 1 & 0 \\ -\cos \chi & 0 & \sin \chi \end{pmatrix} \begin{pmatrix} \cos \psi & -\sin \psi & 0 \\ \sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \end{pmatrix} . \end{aligned} \tag{C.10}$$

In order to align the axis uˆ of the mobile frame with the axis xˆ of the static coordinate system, no rotation about wˆ is needed and thus ψ = 0. The other two rotational angles are defined by the components of the Fermi velocity vector as follows:

$$\eta = \arccos \frac{v\_z}{|v|},\tag{C.11}$$

$$\chi = \begin{cases} \arctan\frac{v\_y}{v\_x}, & v\_x > 0 \\ \arctan\frac{v\_y}{v\_x} + \pi \,, \quad v\_x < 0 \, . \\ \text{sgn}(v\_y)\frac{\pi}{2} \,, \qquad v\_x = 0 \end{cases} . \tag{C.12}$$

Using this transformation, the vector potential and gap parameter can be expressed in the variables of the mobile frame: A(r(X, Y, Z)) and ∆(r(X, Y, Z)). On a characteristic line defined by the impact parameters Y<sup>p</sup> and Zp, the functions from Eq.(C.1) become [257]

$$
\Delta(X) = \Delta(r(X, Y\_p, Z\_p))\,, \tag{\text{C.13}}
$$

$$i\ddot{\epsilon}\_n(X) = i\epsilon\_n + \upsilon\_F eA(r(X, Y\_p, Z\_p))\,,\tag{C.14}$$

$$
\hat{g}(X) = \hat{g}(\mathbf{r}(X, Y\_p, Z\_p), \mathbf{p}\_F, i\epsilon\_n). \tag{C.15}
$$

The parametrisation of the Eilenberger equations on this 1D line is called the *Riccati parametrisation*. Eq. (C.1) formally reduces to the solution of two initial value problems, where the differential equations are scalar and of first order. The Eilenberger propagator is parametrised in terms of two scalar complex functions a(X) and b(X):

$$\hat{g}(X) = \frac{1}{1 + a(X)b(X)} \begin{pmatrix} 1 - a(X)b(X) & 2ia(X) \\ -2ib(X) & -1 + a(X)b(X) \end{pmatrix} . \tag{C.16}$$

The differential equations that need to be numerically solved for a(X) and b(X) are [259]

$$
\hbar v\_F a'(X) + [2\ddot{\epsilon}\_n + \Delta^\dagger(X)a(X)]a(X) - \Delta(X) = 0,\tag{C.17}
$$

$$
\hbar v\_F b'(X) - [2\ddot{\epsilon}\_n + \Delta(X)b(X)]b(X) + \Delta^\dagger(X) = 0 \tag{C.18}
$$

with boundary conditions

$$a(-\infty) = \frac{\Delta(-\infty)}{\epsilon\_n + \sqrt{\epsilon\_n^2 + |\Delta(-\infty)|^2}},\tag{C.19}$$

$$b(+\infty) = \frac{\Delta^\dagger(+\infty)}{\epsilon\_n + \sqrt{\epsilon\_n^2 + |\Delta(+\infty)|^2}}.\tag{C.20}$$

**Figure C.1: Riccati Parametrisation**: A sketch of the calculation procedure using the Riccati parametrization: The 3D trajectories (black arrow) through the flux tube (blue) are given by the Fermi velocities for each point k<sup>F</sup> (θ, φ) on the Fermi surface. Solution of the effective 1D Eilenberger equations for this trajectory yields the LDOS(x) along this line. For a complete picture, a summation over all k<sup>F</sup> (θ, φ) and impact parameters y is performed.

In order to solve them at a certain energy E < ∆0, the analytical continuation iϵ<sup>n</sup> → E + i0 <sup>+</sup> is used. Finally, the local density of states at r is obtained from the real part of g(r(X)), i.e.

$$\mathcal{N}(\mathbf{p}\_F) = \mathcal{N}\_0(\mathbf{p}\_F) \text{Re}\left(\frac{1 - a(X)b(X)}{1 + a(X)b(X)}\right). \tag{C.21}$$

In order to obtain the total local density of states, one still needs an integration over the Fermi surface to calculate all trajectories that are present due to the various Fermi velocity vectors that exist, plus an integration over the impact factors Y in order to account for trajectories that do not traverse the vortex centre. An integration over Z is redundant since the solution in the bulk is the same for every Z. Near the surface, this is strictly not the case anymore but the effect is small and with a large enough variety of Fermi velocity vectors (which we have) these trajectories are not missed. The above described calculation procedure in the Riccati parametrisation is schematically sketched in Fig. C.1.

**Figure C.2: Influence of Vector Potential**: The LDOS obtained from solutions of the 3D Eilenberger equations for a single-flux vortex at energy eU = 0.8 meV without (a) and with vector potential (b), as formulated in Eq. (C.6), shows only quantitative differences. With non-zero vector field, the star arms are still split, yet the CdGM states at this energy are squeezed into a smaller area around the vortex core.

### **C.1 Intermediate Results**

The inclusion of a non-zero vector potential A in the calculations increases the average time needed to solve the Eilenberger equations because the Matsubara frequencies in Eq. (C.14) gain a position dependent term that recquires a coordinate transformation between the two reference frames. Therefore, the simulations shown in the main text are performed without vector potential. It is already visible from Eq. (C.6) that in a calculation where A and ∆ are not solved self-consistently, A only enters the equation as an effective energy term. With a vector potential in the azimuthal direction, like in Eq. (C.6), the scalar product with the Fermi velocity is only expected to yield significant contribution for large impact parameters (for Y<sup>p</sup> = 0, v<sup>F</sup> is perpendicular to A). That means trajectories with increasing impact parameter have large LDOS already for smaller distances from the core than in the field free case. The splitting star arms in the LDOS maps should be squeezed to smaller distances from the core. In fact, this is what is seen in the calculations with vector potential at higher energies, as shown in Fig. C.2. This proves that even though there are quantitative differences to the case without vector potential, in the general characteristics, the LDOS patterns remain unchanged.

# **D Supplementary Pb VortexData**

## **D.1 Bias and Magnetic Field Dependence of CdGM LDOS**

Under reversal of bias voltage sign, i.e. the tunnelling current direction (I ↔ −I), or magnetic field direction (Beˆ<sup>z</sup> ↔ −Beˆz), the differential conductance maps show no qualitative change. As a point of proof, Fig. D.1(a,b) shows a normal vortex stabilized at negative field after exposing the sample to −85 mT for positive voltages and (c,d) an anomalous vortex at 18 mT for negative voltages. This illustrates the electron-hole symmetry of the CdGM states and justifies treating them as excitations of a BCS ground state, i.e. describing their time-evolution by a mean-field BdG Hamiltonian.

**Figure D.1: Reversal of Current/Magnetic Field Direction**: At a magnetic field in −z direction (a,b) or at a negative bias voltage (c,d), the LDOS pattern of normal and anomalous vortices does not change qualitatively compared to the patterns shown in Chap. 6.

## **D.2 Self-consistent Solution of the Gap Equation**

The self-consistent solution of the gap equation for a vortex with isotropic Fermi velocity in Fig. D.2 shows that a Kramer-Pesch effect is to be expected in this case: For a radius from the vortex centre of r < 2ξ<sup>0</sup> the recovery of ∆(r) significantly deviates from the function ∆<sup>0</sup> tanh(r/ξ0). The slope of ∆(r) near the centre is roughly twice as steep as expected from a tanh behaviour with universal ξ0. This leads to a core size ξ (c) = h limr→<sup>0</sup> d∆(r) dr i<sup>−</sup><sup>1</sup> that is roughly half as big.

**Figure D.2: Kramer-Pesch Effect**: The self-consistent calculation of the pair potential ∆ for a vortex with isotropic Fermi velocity exhibits a shrinking of its core size according to the Kramer-Pesch effect, i.e. a steeper recovery of ∆ close to the vortex centre that does not follow tanh(r/ξ0).

# **E Supplementary PbDefectData**

**Figure E.1: SFT in Normal Phase**: dI/dU map (eU = ∆2) of the defect in Fig.7.4(a,e) at B<sup>z</sup> = 100 mT (normal conducting phase). No QPI is visible.

In the normal state at B > Bc, the SFT does not exhibit the characteristic LDOS pattern known from the superconducting state at eU = ∆2, as Fig. E.1 shows. A QPI pattern is also missing.

In the superconducting state, without pinned flux, the SFT misses the zero energy edge state. As Fig. E.2 shows, the superconducting gap is free of states in this region.

**Figure E.2: Missing Edge State**: dI/dU maps at U = 0.0 mV (a) and U = 0.2 mV (b) on the SFT from Fig. 7.4(c) in the superconducting state at B = 31 mT with no pinned vortex. An in-gap edge state of the SFT like in Fig. 7.8(a) and Fig. 7.9(a) is missing.

# **Acknowledgements**

Of course, this thesis would not have been possible without the support, encouragement and guiding by many wonderful people. I am deeply grateful to my supervisor **Prof. Wulf Wulfhekel** for giving me the opportunity to do my doctorate on such an interesting topic. I highly appreciate his continuous support as well as the scientific discussions, the supervision and the shared ideas on such a frequent basis.

I would like to thank **Prof. Jörg Schmalian** for being the second referee, for many helpful discussions and the collaboration on joint publications.

Thanks also go to our staff for always having our backs, so we can concentrate on our research. Special thanks go to **Steffi Baatz**, **Claudia Alaya**, **Michael Meyer** and **Jannis Ret**.

I have to thank my colleague **Dr. Qili Li** for his tireless effort both in the lab and in computational matters. I also thank **Dr. Rolf Heid** for the DFT calculations and many fruitful discussions. Speaking of discussions, I am very grateful to **Dr. Roland Willa**, **Dr. Ioan Pop** and **Prof. Alexander Shnirman** for the discussions and their help in related fields of physics. Special thanks also go to **Dr. Fang Yang** for the hard work on our joint publication and to him and **Dr. Jasmin Jandke** for the great preliminary work.

I want to thank my students **Mirjam Henn** and **Janic Beck** for their interest in my field of study, the great teamwork and the fruitful discussions.

I am lucky to have had many excellent physicists as my teachers throughout the years. I thank you for all your valuable lessons in and out of the lab, **Dr. Timofey Balashov**, **Dr. Toshio Miyamachi**, **Dr. Timo Frauhammer**, **Dr. Marie Hervé** and **Dr. Hung-Hsiang Yang**.

I really enjoyed my stay in the group and want to say thank you to each and everyone that made this time unforgettable. Thank you, **Wulf**, **Timofey**, **Lukas**, **Kevin**, **Marie**, **Fang**, **Sergey**, **Qili**, **Hung-Hsiang**, **Hongyan**, **Loïc**, **Julian**, **Vibhuti**, **Simon**, **Philipp**, **Mirjam**, **Gabriel**, **Qing**, **Haoran**, **Paul** and **Timo**.

Lastly, my sincerest gratitude goes to my wonderful friends and family for their incredible support at each step of the way.

#### **Experimental Condensed Matter Physics (ISSN 2191-9925)**




# Physikalisches Institut Karlsruher Institut für Technologie

Complementary to scattering techniques, scanning tunnelling microscopy provides atomic-scale real space information about a material's electronic state of matter. State-of-the-art designs of a scanning tunnelling microscope (STM) allow measurements at millikelvin temperatures with unprecedented energy resolution. Therefore, this instrument excels in probing the superconducting state at low temperatures and especially its local quasiparticle excitations as well as bosonic degrees of freedom.

Here, three very different types of superconducting materials that are constant subject of fundamental research are studied with a millikelvin STM: Granular aluminium - a granular superconductor, Y123 and Bi2212 - two unconventional cuprate superconductors, and lead (Pb) - a conventional multiband superconductor. In the presence of local inhomegeneities, such as localized spins, vortices or crystal defects, the local quasiparticle excitation spectrum revealed Yu-Shiba-Rusinov states, Caroli-de Gennes-Matricon states and interband coupling effects. Fingerprints of collective excitations, such as antiferromagnetic spin fluctuations or plasmons, were found outside of the superconducting gap and could be used to gain insights into the pair forming but also pair breaking mechanism in the superconductor. Results in this work show that granular aluminium intrinsically hosts localized spins, which raises concern about the performance limit of aluminium/aluminium-oxide based superconducting devices and that (multiple-quanta) vortices are stable in the intermediate state of bulk Pb, which turns out to be an optimal material to study extended vortex core states in combination with multiband superconductivity.

Thomas Gozlinski

**Local Probing of a Superconductor's QPs and Bosonic Excitations with an STM**

**32**

ISSN 2191-9925 ISBN 978-3-7315-1279-0